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ABSTRACT 


Advanced  mechanical  surface  enhancement  techniques  (SET)  have  been  used  successfully  to 
increase  the  fatigue  life  of  metallic  components.  These  techniques  impart  deep  compressive 
residual  stresses  into  the  component  to  counter  potentially  damage-inducing  tensile  stresses 
generated  under  service  loading.  Laser  peening  (LP)  is  an  advanced  mechanical  SET  used 
predominantly  in  the  aircraft  industry.  To  reduce  costs  and  make  the  technique  available  on  a 
large-scale  basis  for  industrial  applications,  simulation  of  the  LP  process  is  required.  Accurate 
simulation  of  the  LP  process  is  a  challenging  task,  because  the  process  has  many  parameters  such 
as  laser  spot  size,  pressure  profile,  and  material  model  that  must  be  precisely  determined.  In  the 
LP  process  material  is  subjected  to  strain  rates  up  to  1  ()6.v  which  is  very  high  compared  to 
conventional  strain  rates.  The  importance  of  an  accurate  material  model  increases  because  the 
material  behaves  significantly  different  at  such  high  strain  rates.  One  of  the  objectives  of  this 
research  is  to  make  advancements  in  the  simulation  of  residual  stresses  induced  by  LP.  Validation 
of  various  material  models  under  investigation  that  could  be  used  in  simulation  and  design  is 
performed.  Inverse  optimization-based  methodology  is  developed  for  simulation  of  residual 
stresses  for  materials  such  as  Inconel®718.  The  procedure  involves  optimizing  the  model 
constants  for  one  load  case  and  using  the  same  constants  for  other  load  cases.  The  second  aspect 
of  this  research  is  to  develop  a  framework  for  uncertainty  quantification  of  the  residual  stress  field 
induced  by  the  LP  process  by  propagation  of  regression  uncertainty.  Development  methodology 
includes  identification  of  regression  uncertainty  as  a  source  of  input  uncertainty  and  using  the 
bootstrap  method  to  verify  the  multivariate  normality  assumption  of  the  model  constant  estimates. 
The  propagation  of  the  input  uncertainty  is  performed  using  Taylor  series  expansion  and 
sensitivity  analysis.  A  confidence  band  for  the  entire  residual  stress  field  is  obtained  and  validated 
using  the  Monte  Carlo  analysis. 
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PREFACE 


This  is  the  final  volume  of  technical  reports  produced  for  the  project  Laser  Peening  for  reliable 
fatigue  life ,  Contract  No.FA8650-04-D-3446-25.  This  volume  covers  the  inverse  optimization  of 
high  strain  rate  properties  from  experimental  material  data,  material  model  validation,  and 
development  of  a  framework  for  uncertainty  quantification  of  residual  stress  fields. 

The  first  volume  included  the  3D  simulation  of  LP,  showing  the  important  process  parameters,  an 
optimization  of  the  residual  stress  field  from  these  key  parameters,  an  outline  of  the  optimization 
strategies  applied,  and  finally  finding  multiple  optima  using  a  modified  particle  swarm 
optimization  method. 
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1.0  SUMMARY 


Laser  Peening  (LP)  has  been  advancing  as  an  able  substitute  for  conventional  treatments  in  the 
process  of  improving  fatigue,  corrosion,  and  wearing  resistance  of  metals.  The  favorable 
compressive  residual  stresses  developed  by  LP  can  extend  further  below  the  surface  than  the 
residual  stresses  from  shot  peening.  Also,  it  is  well  suited  for  precisely  controlled  treatments  of 
localized  fatigue,  and  critical  areas  such  as  holes,  notches,  and  small  fillets  which  might  be 
inaccessible  through  shot  peening.  The  applicability  of  LP  on  complicated  geometries  provides  a 
unique  advantage  when  compared  with  other  surface  enhancement  techniques  (SET)  because 
laser  beams  can  easily  reach  intricate  locations. 

While  the  processing  cost  of  LP  is  high  compared  to  its  counterparts,  advances  in  simulation 
accuracy  and  capabilities  will  decrease  the  simulation  costs  considerably  translating  to  reduction 
in  processing  costs  through  decreased  design  times.  An  ideal  simulation  model  will  incorporate 
the  laser  beam  parameters,  spatial  profiles,  overlay  conditions,  and  material  properties  to  increase 
comparability  between  simulated  and  experimental  results.  The  basic  effect  of  process  parameters 
has  been  extensively  investigated  and  discussed  in  Volume  1.  This  report  details  the  effect  of 
material  behavior  and  the  uncertainties  associated  with  it  when  obtaining  realistic  results  from  the 
simulations. 

Material  models  help  effectively  characterize  the  material  behavior  in  Finite  Element  Analysis 
(FEA).  Three  material  models  have  been  considered  to  validate  the  simulation  results  of 
experimental  results  found  in  literature.  The  models  considered  are  the  Elastic  Perfectly  Plastic 
(EPP),  Johnson-Cook  (JC),  and  Zerilli-Armstrong  (ZA)  material  models.  It  was  found  that  for  the 
various  conditions  the  JC  model  performed  best  of  the  three  models  investigated. 

In  most  cases,  there  is  little  experimental  data  available  to  compare  to  simulated  results.  Using  a 
technique  based  on  inverse  optimization,  the  material  behavior  was  obtained  by  optimizing  the 
material  model  constants  and  reducing  the  error  in  predicting  residual  stresses.  Two  different 
materials  were  tested  to  validate  this  approach.  The  advantage  of  inverse  optimization  is  that 
testing  for  material  model  constants  at  strain  rates  consistent  with  those  experienced  in  the 
simulations  is  that  further  material  model  constant  testing  is  not  necessary. 

This  research  also  develops  a  framework  to  determine  the  uncertainty  in  the  resultant  residual 
stresses  from  the  LP  simulation.  Material  model  constants  are  considered  to  be  uncertain.  The 
uncertainty  represents  the  variability  of  the  material  model  when  characterizing  the  material 
behavior.  The  ’Bootstrapping  Method’  employed  in  statistics  is  used  for  this  process.  Two  case 
studies  were  used  to  validate  this  method.  Finally,  the  method  was  applied  to  an  LP  simulation, 
cementing  the  suitability  of  the  JC  model  for  LP  simulations. 
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2.0  INTRODUCTION 


This  chapter  introduces  the  fatigue  life  issues  of  the  aircraft  lug  component  and  the  motivation  for 
the  development  of  SET.  An  introduction  to  SET  is  presented.  The  mechanism  of  the  SET  is 
described,  which  generates  favorable  compressive  residual  stresses  in  the  surface  regions.  Several 
techniques  such  as  shot  peening,  Low  Plasticity  Burnishing  (LPB),  waterjet  peening,  and  LP  are 
discussed. 

2.1  Motivation 

The  F-22  Raptor  is  a  fighter  aircraft  that  uses  stealth  technology.  It  is  primarily  an  air  superiority 
fighter.  The  United  States  Air  Force  considers  the  F-22  a  critical  component  of  the  fleet  [1], 

Figure  1  shows  a  picture  of  two  F-22  raptors  in  column  flight  formation  [2]  . 


Figure  1:  Two  F-22  Raptors  in  Column  Flight  Formation  [2] 


To  demonstrate  the  capability  of  the  design,  the  F-22  aircraft  design  was  subjected  to  a  full-scale 
fatigue  test  by  Cayton  and  his  co-workers  [3].  The  main  goal  of  this  full-scale  testing  was  to 
identify  any  anomalies  that  remained  undiscovered.  One  such  structural  problem  was  found  in  the 
form  of  a  crack  on  one  of  the  lugs  of  the  wing  attachments.  Figure  2  shows  the  location  of  the 
crack  on  the  aircraft  [4].  The  lower  lugs  are  important  for  durability  and  damage  tolerance  [3]. 
The  summary  of  the  analysis  was  that  the  crack  was  a  result  of  ForceMate™  expansion  levels  for 
the  applied  load  levels.  The  detailed  report  is  available  in  the  work  of  Cayton  et  al.  [3].  A 
configuration  change  was  made  in  the  design  of  the  aircraft  and  the  change  was  implemented  in 
the  manufacturing  process  for  subsequent  aircraft.  Since  there  were  aircraft  already  manufactured 
with  the  old  design,  a  program  has  been  proposed  to  increase  the  fatigue  life.  This  program 
suggests  to  use  techniques  to  increase  the  fatigue  life  without  replacing  the  component.  SET  are 
one  possible  option  that  can  address  the  issue. 
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2nd  LP  Conference,  2010 


Figure  2:  Location  of  Crack  on  Lug  of  an  Aircraft  Lug  [4] 


2.2  SET 

Highly  stressed  components  such  as  aircraft  lugs,  turbine  blades  from  aerospace  applications, 
connecting  rods  from  automobiles,  and  knee  implants  from  medical  applications  demand  an 
optimal  judgment  of  material,  loading  conditions,  design,  and  production  processes  in  the 
conceptual  stages  (Figure  3). 


Conceptual  Stage 


Figure  3:  Attributes  Considered  During  Conceptual  Stage  of  a  Component 
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The  surface  of  the  component  plays  an  important  role  under  fatigue  loading  because  cracks  are 
initiated  from  the  surface  discontinuities  such  as  notches,  holes,  and  voids,  ultimately  leading  to 
failure.  With  respect  to  this,  various  SET  are  being  used  in  the  industry  to  mitigate  or  delay 
failure.  These  techniques  can  be  broadly  classified  [5]  as: 

1.  Mechanical  surface  treatments 

2.  Surface  diffusion  treatments 

3.  Surface  overlay  coatings 

The  basic  idea  of  mechanical  treatments  is  to  induce  compressive  residual  stresses  on  the  surface 
regions.  Residual  stresses  are  defined  as  stresses  which  exist  in  the  bulk  of  the  material  without 
application  of  external  load.  These  compressive  stresses  negate  the  tensile  stresses  acting  on  the 
surface  leading  to  increased  life  of  the  component.  Figure  4  shows  the  basic  mechanism  for  the 
production  of  residual  stresses  in  mechanical  SET. 


After  SET 
L+AL 

< - ► 

— ►  A  <•— 


Figure  4:  Basic  Mechanism  to  Generate  Residual  Stresses 


Before  any  pressure  is  applied  to  the  surface,  the  length  of  surface  and  core  are  the  same.  When  a 
high  pressure  is  applied  on  the  surface,  shock  waves  are  generated  and  progress  through  the 
depth.  This  causes  plastic  deformation  in  the  top  region.  Since  the  magnitude  of  the  shock  wave 
decreases  as  it  moves  down  to  the  core,  it  does  not  produce  plastic  deformation  in  the  core.  To 
keep  cohesion  between  the  stretched  surface  and  the  core,  the  layers  in  the  surface  are  set  to 
compressive  stresses.  To  maintain  equilibrium,  compensatory  tensile  stresses  are  induced  in  the 
core.  In  this  section  several  of  the  techniques  used  in  the  industry,  such  as  shot  peening,  LP,  LPB, 
and  waterjet  peening  are  discussed. 

2.2.1  Shot  Peening 

Shot  peening  is  one  of  the  traditional  techniques  used  to  induce  compressive  residual  stresses. 
The  discovery  of  improvement  in  fatigue  properties  is  seen  from  efforts  in  1928  and  1929  by 
Buick  Motor  Division  [6].  Shot  peening  is  a  cold  working  process  in  which  the  surface  of  the 
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component  is  bombarded  with  small  spherical  media  called  shot  [7].  Each  piece  of  shot  striking 
the  surface  acts  as  a  tiny  peening  hammer,  imparting  to  the  surface  a  small  indentation  or  dimple. 
A  schematic  of  the  process  is  shown  in  figure  5. 


Shot 


Dimple 

A 


V- 


d 


Figure  5:  Schematic  of  the  Shot  Peening  Process 


Steel  balls  are  the  popular  material  used  for  the  shot.  The  most  commonly  used  ball  size  is  0.05  to 
1  mm  in  diameter.  Shot  peening  has  the  advantage  of  low  cost  and  has  a  successful  track  record  of 
application  to  compressor  blades  [8].  A  three  dimensional  finite  element  simulation  was 
performed  by  Meguid  et  al.  using  ANSYS  [9].  Three  main  issues  addressed  in  this  work  were:  (i) 
effect  of  shot  velocity,  size,  and  shape  upon  the  plastic  zone  development,  (ii)  effect  of  the 
separation  distance  between  two  shots  on  the  unloading  residual  stresses,  and  (iii)  effect  of 
strain-hardening  rate  of  the  target  upon  the  spread  of  the  plastic  zone.  Shot  peening  has  been  used 
in  the  industry  during  the  past  decade.  The  other  advantage  of  using  shot  peening  is  that  the  cost 
involved  is  less  compared  to  other  techniques.  The  process  also  has  disadvantages.  The  surface 
finish  is  damaged  due  to  roughening  by  the  balls.  The  depth  of  compressive  residual  stresses 
generated  by  shot  peening  is  relatively  low  compared  to  other  techniques.  This  is  considered  to  be 
one  of  the  major  disadvantages  of  the  method  [10]. 

2.2.2  LPB 

LPB  is  a  SET  developed  to  induce  compressive  stresses  with  higher  depth  compared  to  shot 
peening  and  minimal  cold  work  [11].  In  this  process,  a  smooth,  free  rolling  ball  makes  a  single 
pass  over  the  material  under  a  normal  force  just  sufficient  to  plastically  deform  it.  This  creates 
compressive  residual  stresses  on  the  surface  region.  The  process  is  shown  schematically  in  Figure 
6. 
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Normal  Force 


Figure  6:  Schematic  of  the  LPB  Process 


The  ball  is  supported  in  a  fluid  bearing  with  sufficient  pressure  to  lift  the  ball  off  of  the  surface  of 
the  retaining  spherical  socket.  The  ball  is  in  solid  contact  only  when  the  surface  needs  to  be 
burnished.  It  is  free  to  roll  otherwise.  The  force  is  applied  from  top  and  the  tool  position  is 
generally  computer  controlled  [12].  LPB  has  been  used  in  medical  field  for  hip  stems  [13].  The 
LPB  process  was  developed  and  applied  to  the  modular  neck  taper  junction  of  total  hip  prosthesis. 
The  advantages  of  LPB  are  its  low  cost  and  minimal  cold  work.  The  disadvantage  of  this 
technique  is  that  it  cannot  be  applied  to  complicated  geometries  such  as  a  mechanical  gear. 

2.2.3  Waterjet  Peening 

Waterjet  peening  is  very  similar  to  shot  peening,  but  instead  of  using  balls,  water  with  high 
velocity  is  impinged  upon  the  target  [14].  A  schematic  of  the  process  is  shown  in  Figure  7. 
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Figure  7:  Schematic  of  the  Water  jet  Peening  Process 


The  waterjet  generates  a  pressure  distribution,  creating  localized  plastic  deformation.  This 
deformation  is  restrained  by  surrounding  material,  leading  to  compressive  residual  stress  in  the 
surface  regions.  The  advantages  of  waterjet  peening  compared  to  shot  peening  are  a  better  surface 
finish  and  increased  coverage  capabilities.  The  disadvantage  is  that  the  magnitude  of  residual 
stress  is  lower  than  other  treatments  and  low  controllability.  There  are  various  parameters  such  as 
water  pressure,  nozzle  feed  speed,  peening  duration,  nozzle  and  peening  angle  that  control  the 
magnitude  of  residual  stresses  [15].  Various  approaches  [16,  17]  are  available  for  mathematical 
modeling  of  the  process  that  vary  from  closed  form  differential  equations  to  finite  element 
method. 

2.2.4  LP 

LP  is  a  very  high  powered  process  with  a  short  interactive  time.  It  requires  a  laser  power  density 
on  the  order  of  109!F cm~2  with  pulse  duration  in  nano  seconds.  Figure  8  shows  a  schematic  of 
the  process. 
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Shock  Waves 


^^Transparent  Overlay 


Figure  8:  Schematic  of  the  LP  Process 


In  the  LP  ,  a  laser  beam  is  fired  at  an  opaque  overlay  (example:  black  paint)  applied  to  the  surface 
of  the  component  to  be  laser  peened.  The  opaque  overlay  is  covered  in  turn  with  a  transparent 
overlay  such  as  water.  This  water  overlay  can  either  be  running  or  still  water.  The  vaporization  of 
the  opaque  overlay  produces  a  plasma  confined  between  the  two  overlays.  The  pressure  exerted 
on  the  material  surface  by  the  expanding,  confined  plasma  generates  shock  waves  that  propagate 
into  the  target.  Depending  on  the  laser  parameters,  if  the  magnitude  of  shock  waves  is  high 
enough,  the  material  plastically  yields,  resulting  in  compressive  residual  stresses  in  the  surface 
regions.  The  formation  of  plastic  strain  by  shock  waves  continues  until  the  peak  stress  wave 
decreases  below  the  dynamic  yield  strength.  Compensatory  tensile  stresses  are  also  created  inside 
the  target.  Details  about  the  process  can  be  found  in  the  review  papers  of  Peyre  and  Fabbro  [18] 
and  Mantross  et  al.  [19].  The  advantage  of  the  LP  process  is  that  it  can  be  performed  on  complex 
geometries  and  involves  no  contact  with  the  component.  Compared  to  the  shot  peening  process, 
LP  achieves  higher  depth  of  compressive  residual  stress. 

2.3  Section  Summary 

Fatigue  life  issues  in  engineering  structures  are  introduced  with  focus  on  aircraft  component.  A 
basic  mechanism  of  SET  is  introduced.  Various  SET  are  briefly  described  including  advantages 
and  disadvantages. 
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3.0  BACKGROUND 


This  chapter  provides  the  background  details  on  the  LP  process.  Evolution  of  the  LP  process  in  a 
historical  perspective  is  described.  A  review  of  LP  experiments  performed  on  different  materials 
is  provided.  The  parameters  involved  in  the  LP  process  for  a  successful  implementation  are 
discussed.  A  review  of  the  simulation  of  the  LP  process  is  discussed  and  the  advantages  and 
limitations  of  the  work  are  provided.  Also,  a  review  of  uncertainty  quantification  methodologies 
relevant  to  spatial  structural  response  such  as  residual  stresses  is  detailed.  The  limitations  of  the 
review  leads  to  the  definition  of  research  need. 

3.1  Historical  Perspective 

The  history  and  evolution  of  the  LP  process  has  been  described  in  detail  by  Clauer  [20].  Clauer 
divided  the  time  between  the  discovery  of  the  phenomenon  in  1960,  to  the  phase  into  production 
in  the  1990  into  stages.  The  first  experimental  generation  of  shock  waves  in  a  laboratory  was 
performed  by  White  [21]  in  1962.  The  shock  waves  were  generated  by  electron  bombardment  and 
electromagnetic  wave  absorption.  The  experiments  were  conducted  in  a  vacuum  chamber  and 
peak  power  densities  of  44 W / cm 2  were  achieved.  Gregg  and  Thomas  [22]  and  Skeen  and  York 
[23]  were  a  few  of  the  researchers  who  extended  the  investigations  on  characterizing  the 
momentum  transfer  phenomena  for  the  laser  induced  plasma.  The  enhanced  pressure  was 
achieved  by  the  use  of  the  transparent  overlay  was  discovered  by  Anderholm  in  the  1970  [24]. 
Stress  measurements  were  made  using  quartz  gauze.  Aluminum  foil  and  quartz  disk  were  used  as 
overlays.  Although  the  experiment  was  performed  in  a  vacuum  chamber,  the  pressure 
enhancement  can  also  be  achieved  in  air.  The  benefits  of  laser- induced  shock  waves  in  material 
properties  were  identified  by  Malozzi  and  Lairand,  who  were  awarded  the  patent  in  1974  [25]. 
There  were  additional  patents  awarded  until  the  present  date  in  various  applications  of  LP. 

During  the  1970  and  early  1980  efforts  were  made  to  model  laser  material  interaction  using  LILA 
and  TOODY  codes.  Many  of  the  efforts  were  made  towards  understanding  the  effects  of  laser 
processing  parameters  on  fretting  fatigue,  hardening,  galling,  and  wear  of  metals.  Clauer  and 
Lairand  demonstrated  that  compressive  residual  stresses  were  a  contributing  factor  to  the  increase 
in  fatigue  life  of  metallic  components  [26] .  The  first  step  into  production  was  made  in  early 
1990’s,  when  LP  was  used  on  several  L101  fan  blades  to  reduce  the  foreign  object  damage.  The 
goal  at  the  current  stage  is  to  lower  the  processing  cost  and  to  engineer  the  process  parameters  for 
new  applications. 

3.2  Experimental  Attribute  of  LP 

The  LP  process  has  three  major  ingredients  for  a  successful  application.  They  are  (i)  laser  system, 
(ii)  opaque  overlay,  and  (iii)  transparent  overlay.  This  section  gives  a  brief  description  about  each 
of  the  three. 

3.2.1  Laser  Systems 

The  laser  system  for  the  LP  process  is  not  one  of  the  standard  laser  systems  available  in  the 
industry.  These  system  have  to  be  custom  built  because  the  minimum  energy  required  is  20 


9 

Approved  for  public  release;  distribution  unlimited. 


Joules/pulse,  whereas  the  standard  in  the  industry  is  3  Joules/pulse.  The  wavelength  of  the  system 
is  very  important  because  it  controls  the  interaction  with  the  target  to  generate  the  shock  waves. 
The  wavelength  should  be  about  lj am.  These  constraints  reduce  the  large  variety  of  laser  systems 
available  to  two  or  three  viable  possibilities.  Pulsed  CO2  laser  systems  are  proven  to  be  efficient 
and  durable,  but  they  have  a  characteristic  wavelength  of  10.6 /dm  that  makes  them  infeasible  for 
the  LP  process.  Currently,  neodymium-glass  laser  is  the  most  suitable  candidate.  Detailed 
description  of  the  selection  of  laser  systems  is  provided  by  Shepard  [27]  and  Fairand  and  Clauer 
[28].  The  experiments  performed  by  Fairand  and  Clauer  [29]  used  a  Q-switched 
neodymium-glass  laser  that  consisted  of  an  oscillator  followed  by  six  amplifier  stages.  The 
system  was  capable  of  emitting  up  to  500  J  of  laser  energy  in  a  pulse  with  a  full  width  at  half 
maximum  (FWHM)  of  20-40  ns.  Peyre  et  al.  [30]  used  a  Nd-glass  prototype  laser  operating  at  a 
wavelength  of  1.06 pm  followed  by  four  amplifier  stages.  This  system  is  capable  of  emitting  up  to 
80  J  in  a  pulse  that  is  semi-gaussian  and  short  rise  time  (SRT)  in  shape  with  a  FWHM  of  15-30  ns. 

3.2.2  Opaque  Overlay 

The  opaque  overlay  in  the  LP  process  serves  two  purposes.  First,  vaporization  of  the  opaque 
overlay  forms  plasma  that  generates  the  shock  waves  inside  the  target,  creating  favorable 
compressive  residual  stresses.  Second,  it  acts  as  a  protective  layer  for  the  target.  Various  materials 
such  as  Al,  Zn,Pb,  and  black  paint  have  been  used  as  opaque  overlays.  It  has  been  shown  by 
Fairand  and  Clauer  [29]  that  all  opaque  overlays  produce  nearly  the  same  shock  wave  pressures 
after  crossing  certain  laser  power  density.  The  thickness  of  the  black  paint  is  typically  8  —  10 pm. 
Self-adhesive  foils  were  proven  to  be  better  in  certain  conditions  compared  to  protective  coatings 
[31]  The  practical  advantage  of  black  paint  compared  to  others  has  made  it  attractive  in  the 
commercial  LP  process. 

3.2.3  Transparent  Overlay 

The  presence  of  transparent  overlay  is  known  as  confined  regime,  compared  to  direct  regime, 
which  does  not  have  any  confining  overlay.  The  earlier  experiments  that  induced  shock  waves  by 
laser  systems  did  not  include  transparent  overlays.  Transparent  overlay  was  first  used  by 
Anderholm  [32,  24].  The  intensity  of  shock  waves  due  to  confinement  by  transparent  overlays  is 
up  to  two  orders  of  magnitude  higher  than  plasmas  generated  in  vacuum  [33,  34,  35].  According 
to  Peyre  et  al.  [36],  the  pressure  of  shock  waves  is  five  to  ten  times  higher  in  magnitude  and  two 
to  three  times  longer  in  duration  than  direct  regime.  Fused  quartz,  glass,  polymeric  material  and 
distilled  water  have  been  used  successfully  as  a  transparent  overlay.  Water  is  the  most  widely 
used  in  the  commercial  LP  process.  This  water  can  be  either  running  water  or  stagnant  water. 

3.3  Review  of  Simulation  of  the  LP  Process 

One  of  the  major  disadvantage  of  the  LP  process  is  the  cost  involved.  The  fixed  cost  to  set  up  the 
laser  system  is  high,  and  the  variable  operating  costs  for  particular  applications  are  also  high.  One 
reason  for  higher  operating  costs  can  be  attributed  to  the  current  industry  practice  of  performing 
LP  experiments  and  then  obtaining  the  optimum  process  parameters  for  each  application. 
Simulation  of  the  LP  process  not  only  reduces  the  requirements  for  expensive  trial-and-error 
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experiments,  it  also  helps  in  engineering  the  process  parameters.  The  review  of  the  simulation  of 
the  LP  process  is  described  here  in  chronological  order. 

The  first  finite  element  simulation  to  predict  residual  stresses  from  the  LP  process  was  reported  in 
1999  by  Braisted  and  Brockman  [37].  The  methodology  of  sequential  switch  from  explicit  to 
implicit  algorithms  for  efficient  computation  was  developed.  They  also  introduced  finer  and 
infinite  elements  to  take  advantage  of  the  localized  effect  of  the  LP  process.  Constitutive  behavior 
of  elastic  perfectly  plastic  material  behavior  was  assumed.  Nam  and  co-workers  from  The  Ohio 
State  University  developed  their  own  finite  element  code  called  SHOCKWAVE  to  model  the 
residual  stresses  induced  by  the  LP  process  [38].  The  JC  material  model  was  used  for  the 
constitutive  behavior.  Single-side  and  double-side  LP  were  modeled.  The  shock  wave  profile  was 
measured  by  using  streak  cameras  and  quartz  gages  by  LSP  Technologies,  LLC.  The  shock  wave 
has  temporal  and  spatial  profile.  Experimental  results  were  provided  by  LSP  Technologies. 
Research  efforts  to  model  the  LP  phenomenon  were  carried  out  independently  in  Saudi  Arabia, 
France,  and  Australia  in  2003.  Researchers  in  Saudi  Arabia  worked  on  two  aspects  of  modeling 
the  LP  process.  The  first  one  included  developing  analytical  equations  for  temperature  during 
plasma  generation,  recoil  pressure  calculations,  and  wave  analysis  for  plastic  deformation  [39]. 
The  second  work  focused  on  numerical  prediction  of  the  depth  of  plastic  deformation  and 
compressive  residual  stresses  [40].  Finite  difference  scheme  was  used  to  simulate  the  propagation 
of  stress  waves  coupled  with  finite  element  models.  Peyre  et.  al  from  France  used  ABAQUS  for 
the  prediction  of  residual  stresses.  The  JC  model  with  Equation  of  State  (EOS)  was  used  as  the 
procedure  for  material  behavior.  Axisymmetric  mesh  with  infinite  elements  was  used  to  represent 
the  structure.  Velocity  Interferometer  System  for  Any  Reflector  technique  was  used  to  predict  the 
loading  conditions.  Ding  and  Ye  developed  the  first  3D  simulations  for  prediction  of  residual 
stresses  [41].  They  also  developed  the  two-sided  LP  process  [42].  EPP  constitutive  model  was 
assumed  for  the  material  behavior.  History  of  energies  was  used  as  a  stopping  criteria  for 
ABAQUS/Explicit.  The  3D  simulations  were  also  modeled  by  Ocana  and  his  co-workers  [43]. 
Multiple  overlapping  shots,  which  were  neglected  in  the  previous  work  were  modeled  using 
ABAQUS.  Wu  and  Shin  from  Purdue  University  developed  a  complete  model  for  simulating  the 
residual  stresses  induced  by  the  LP  process.  Their  first  work  was  published  in  2004  [44],  and  their 
work  continued  until  recently  in  2007  [45].  Most  of  the  work  involved  developing  analytical 
equations  for  converting  the  laser  pulse  to  pressure  pulse.  According  to  them  the  plasma 
generated  from  the  LP  process  is  divided  into  breakdown  and  confined  plasma.  Only  the  confined 
plasma  is  able  to  produce  the  shock  waves.  There  is  no  need  for  any  measurement  of  shock  wave 
pressure  magnitudes.  Inputs  are  needed  for  the  laser  equipment  such  as  power  density,  laser 
wavelength.  The  modeling  procedure  predicts  the  residual  stresses.  The  JC  model  was  used  as 
constitutive  model.  The  recent  work  in  modeling  the  residual  stresses  include  Peyre  et  al.  and 
Warren  et  al.  in  2007  [46]  and  2008  [47]  respectively.  The  main  goal  of  Peyre  and  his  co-workers 
was  to  investigate  the  influence  of  protective  coatings,  and  precise  simulation  of 
thermo-mechanical  uncoated  LP  process.  Shock  propagation  and  attenuation  was  modeled  in 
ABAQUS/Explicit,  and  thermal  modeling  of  plasma  heating  was  modeled  in  ABAQUS/Standard. 
Thermo-mechanical  residual  stresses  were  predicted  by  combining  the  results  of 
ABAQUS/Explicit  and  ABAQUS/Standard.  Warren  et  al.  performed  3D  simulations  of  single  and 
multiple  shots.  A  user  subroutine,  VDLOAD,  was  developed  in  ABAQUS  to  include  spatial  and 
temporal  variation.  Spot  size  effect,  spacing  effect,  and  intensity  effect  on  residual  stresses  were 
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investigated.  Srinivasan  [48]  developed  simulation  work  for  the  spallation  response,  caused  by  the 
LP  process  through  strain  rate  and  temperature  dependent  material  models.  The  JC  and  modified 
ZA  models  were  used  as  material  models.  The  constants  were  obtained  from  other  references. 

The  most  recent  work  [49]  involves  3D  simulations  on  curved  surfaces  for  geometric  effects  on 
residual  stress  predictions.  This  is  the  first  effort  to  simulate  a  three  dimensional  curved  geometry. 
There  were  no  experimental  comparisons  in  their  work. 

3.4  Review  of  Relevant  Uncertainty  Quantification  Research 

Uncertainty  quantification  provides  an  alternative  to  conventional  deterministic  analysis.  The 
uncertainty  is  taken  into  consideration  in  deterministic  analysis  by  employing  the  safety  factor 
approach.  This  safety  factor  accounts  for  the  uncertainties  in  the  design,  manufacturing,  and  other 
factors  in  making  the  final  product.  While  deterministic  approach  with  safety  factor  have  been 
used  with  success  in  the  past,  they  tend  to  be  over  conservative.  Today’s  increased  demand  for 
higher  efficiency,  better  performance  and  lower  cost  makes  it  necessary  to  consider  developing 
techniques  to  quantify  these  uncertainties.  Probability  theory  has  been  widely  used  to  quantify  the 
response  caused  by  uncertainties  in  inputs.  Non-probability  based  methods  exist,  such  as  fuzzy 
set  theory  [50,  51],  evidence  theory  [52],  and  information- gap  theory  [53].  Uncertainties  can  be 
broadly  categorized  as  model  uncertainty  and  parametric  uncertainty.  Model  uncertainty  is  caused 
when  there  are  different  models  to  represent  the  same  physical  situation.  Choosing  one  model 
over  other  is  refered  as  model  uncertainty.  Parametric  uncertainty  is  caused  within  a  model.  It  is 
caused  due  to  variations  in  the  input  parameters  for  a  simulation  model.  The  uncertainties  can 
also  be  classified  into  two  categories:  aleatory  and  epistemic.  Aleatory  uncertainty  is  irreducible 
or  inherent  uncertainty.  Epistemic  uncertainty  is  a  result  of  lack  of  knowledge  or  data,  and 
therefore  is  a  reducible  uncertainty.  There  are  many  ways  to  characterize  the  uncertainty  in  the 
inputs  and  propagation  of  the  input  uncertainty  through  a  physics-based  simulation  or  analytical 
equation.  The  final  goal  is  to  predict  the  uncertainties  in  the  response  of  the  systems  caused  by  the 
uncertainties  in  the  inputs. 

Sampling  methods  are  the  most  common  method  to  quantify  uncertainty  in  system  response.  New 
methods  that  are  developed  are  often  compared  with  a  sampling  method  for  verification.  The 
computational  cost  of  sampling  methods  makes  it  impractical  for  implementation.  Monte  Carlo 
Simulation  (MCS)  technique  with  Latin  Hypercube  Sampling  (LHS)  is  a  popular  sampling 
method  [54].  The  drawback  of  the  LHS  is  its  inefficiency  to  sample  for  correlated  input  random 
variables.  This  limitation  has  been  recently  addressed  and  an  extension  to  the  existing  LHS 
technique  is  developed  to  incorporate  correlation  between  the  variables  [55]. 

Techniques  based  on  limit  state  function  were  developed.  These  include  Lirst  Order  Reliability 
Method  (LORM)  and  Second  Order  Reliability  Method  (SORM).  In  LORM,  the  limit  state  is 
approximated  by  linear  approximation  at  the  Most  Probable  Point  (MPP).  The  accuracy  of  the 
LORM  depends  on  failure  surface.  LORM  is  most  accurate  when  the  limit  state  is  nearly  linear  in 
the  vicinity  of  MPP.  LORM  predicts  inaccurate  results  when  the  failure  surface  is  nonlinear.  If  the 
failure  surface  is  approached  by  second  order  approximation  at  the  MPP  it  is  SORM.  The 
addition  of  second  order  term  increases  the  efficiency,  but  also  increases  the  computational  cost. 
The  MPP-based  methods  can  be  based  on  high  quality  function  approximations  [56].  This  is 
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based  on  standard  normalized  system.  Performance  measure  approach  [57]  is  developed  that  is 
shown  to  be  more  stable  and  efficient  than  reliability  index. 

Approximate  integration  methods  have  been  developed  to  estimate  the  statistical  moments  of  the 
response.  Univariate  dimension  reduction  method  was  developed  by  Rahman  and  Xu  [58]  for 
numerical  integration  of  multi-dimensional  integration.  It  converts  a  multidimensional  integration 
into  multiple  one  dimension  integration.  Recently  this  method  has  been  extended  by  Youn  et  al. 
by  including  eigenvector  sampling  technique  among  other  developments  [59].  Penmetsa  and 
Grandhi  have  adapted  the  fast  fourier  transformation  (FFT)  technique  to  estimate  the  structural 
failure  probability  [60].  Convolution  integral  is  converted  into  multiplication  by  use  of  FFT.  Two 
point  adaptive  nonlinear  approximation  is  used  at  the  MPP 

Spatial  variability  is  another  form  of  uncertainty  that  quantifies  variability  in  the  structural 
response.  Choi  et  al  developed  framework  by  considering  random  field  variations  [61].  The  effect 
of  spatial  variation  in  Young’s  modulus  on  the  responses  such  as  displacement  is  shown  in  the 
work.  The  framework  developed  involved  use  of  Karhunen-Loeve  transformation  for  reducing  the 
number  of  input  random  variables.  Polynomial  chaos  expansion  is  used  to  construct 
approximation.  Analysis  of  variance  is  performed  for  significance  test.  Probabilistic  analysis  of 
residual  stress  (spatial  structural  response)  is  performed  by  Millwater  and  his  co-workers  [62,  63]. 
In  the  first  work,  the  residual  stress  is  modeled  as  probabilistic  in  nature  by  a  random  scaling 
factor.  The  main  objective  was  to  identify  the  effect  of  residual  stress  on  probability  of  fracture 
for  a  compressor  disk.  Other  random  variables  that  were  considered  in  the  design  process  were 
initial  crack  size  distribution,  crack  propagation  scatter,  and  stress  scatter.  An  extension  to  the 
work  is  performed  by  assuming  random  residual  stress  field.  Two  analytical  models  were  fit  to  the 
experimental  data  using  nonlinear  regression  analysis.  The  randomness  in  the  estimates  obtained 
were  utilized  further  to  obtain  confidence  bounds  on  the  residual  stress  field.  Khaled  and  Noor 
[64]  examined  the  effect  of  uncertainty  in  the  material  properties  of  steel  associated  with 
martensitic  transformations  on  the  residual  stress  field  induced  by  the  welding  process.  They  used 
a  fuzzy-set  approach  to  quantify  the  uncertainty  in  material  properties.  The  main  conclusion  was 
that  the  variation  in  material  properties  has  a  significant  effect  on  welding  residual  stress  field. 
Grujiccic  et.  al.  [65]  quantified  the  effect  of  material  and  processing  parameters  on  the  magnitude 
and  distribution  of  the  axial  residual  stress  on  gun  barrels.  Advanced  mean  value  method,  which 
is  based  on  limit  state  concept,  was  adapted  to  determine  the  cumulative  distribution  of  residual 
stress.  Sobczyk  and  Trebicki  [66]  proposed  one  to  three  parameters  to  model  the  random  residual 
stress  on  fatigue  crack  growth.  The  probabilistic  model  of  residual  stress  was  included  in  their 
evaluation  of  stress  intensity  factor  and  calculations  of  crack  growth  life.  The  Forman  fatigue 
crack  growth  equation  was  used  to  obtain  the  cycles  to  failure.  Park  and  his  co-workers  [67]  have 
developed  an  approach  based  on  bayesian  statistics  to  quantify  the  model  uncertainty  by  using  the 
measured  differences  between  the  experimental  data  and  simulation  predictions.  Model 
probability,  that  is  defined  as  the  degree  of  belief  that  a  model  is  the  best  approximating  model 
among  a  model  set  is  estimated  using  the  experimental  data.  The  methodology  is  applied  to 
obtain  a  confidence  band  on  the  residual  stress  field  induced  by  the  simulation  of  the  LP  process. 

3.4.1  Sources  of  Uncertainty  in  LP  Simulation 

Different  sources  of  uncertainty  exist  during  the  simulation  of  the  LP  residual  stresses.  Figure  9 
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shows  an  overview  of  these  different  sources.  The  uncertainties  can  be  categorized  as: 

1.  Model  Uncertainty 

2.  Parametric  Uncertainty 
3. Shape  Uncertainty 

4. Regression  Uncertainty 


Figure  9:  Sources  of  Uncertainty  During  LP  Simulation 


Model  uncertainty  is  defined  as  the  uncertainty  involved  in  selecting  the  best  model  from  the 
possible  models  to  represent  the  physical  scenario.  Details  regarding  model  uncertainty  can  be 
found  in  the  work  of  Park  [67].  An  example  of  model  uncertainty  would  be  selecting  the  best 
material  model  from  among  the  JC  model,  ZA  model,  and  EPP  model  to  predict  the  material 
behavior  based  on  limited  experimental  data. 

Parametric  uncertainty  arises  due  to  uncertainty  in  the  input  process  parameters.  Parametric 
uncertainty  can  be  categorized  as  aleatory  or  epistemic,  based  on  the  availability  of  experimental 


14 

Approved  for  public  release;  distribution  unlimited. 


data.  Examples  of  parametric  uncertainty  are:  i)  peak  pressure  of  the  shock  wave,  ii)  laser  spot 
radius,  iii)  geometric  dimensions,  and  iv)  material  properties. 

Shape  uncertainty  is  defined  as  the  uncertainty  in  shape  of  the  input  profiles  that  are  essential  for 
the  LP  simulation.  Examples  of  shape  uncertainty  include  temporal  and  spatial  profiles  of  the 
shock  wave.  The  shape  of  the  laser  spot  acts  as  a  shape  uncertainty  for  3D  models. 

Regression  uncertainty  is  defined  as  the  uncertainty  in  the  model  constant  estimates  obtained 
from  the  fit  of  the  constitutive  material  model  to  the  experimental  data.  Regression  uncertainty 
represents  the  uncertainty  in  the  material  model  used  to  predict  the  material  behavior. 

3.5  Definition  of  Research  Need 

The  review  of  the  research  in  simulation  of  the  LP  process  and  uncertainty  quantification  leads  to 
three  major  research  needs.  These  three  research  needs  are  separated  into  two  groups.  The  first 
group  is  advancements  in  finite  element  simulation  of  residual  stresses  induced  by  the  LP  process, 
and  the  second  group  is  development  of  uncertainty  quantification  framework  for  spatial 
structural  response  such  as  residual  stress. 

Due  to  the  complexity  of  the  process,  development  of  reliable  and  accurate  simulation  that  can  be 
used  in  process  design  is  a  challenging  task.  There  have  been  disagreements  between 
experimental  and  simulated  results  in  previous  work.  The  LP  process  has  many  parameters  that 
need  to  be  precisely  determined  for  an  accurate  finite  element  simulation.  Among  them  are 
pressure  profile,  laser  spot,  and  material  model.  The  first  two  parameters  depend  on  the  laser 
system  being  used  in  experiments.  During  the  LP  process  the  component  experiences  strain  rates 
of  lO6/^.  This  work  emphasizes  the  need  for  an  accurate  material  model  because  the  material 
behaves  significantly  different  at  such  high  strain  rates.  There  has  been  no  work  reported  on 
validation  of  material  model  for  the  simulation  of  the  LP  process. 

Secondly,  most  of  the  materials  used  for  simulation  in  the  previous  work  had  little  experimental 
data  of  the  material  behavior.  This  insufficient  or  little  experiment  data  was  used  to  predict  the 
model  constants.  The  material  model  was  extrapolated  for  higher  strain  rates,  where  there  was  no 
experimental  data  available.  In  this  work,  LP  experiments  on  Inconel®718  were  performed  in 
collaboration  with  LSP  Technologies.  Lor  Inconel®718  material,  no  experimental  data  of  the 
material  behavior  for  room  temperatures  at  high  strain  rates  were  found.  The  simulation  of  the  LP 
process  for  materials  which  don’t  have  any  experimental  data  is  non  trivial. 

Thirdly,  uncertainty  quantification  of  spatial  structural  response,  such  as  residual  stress  from  the 
LP  process  is  needed.  No  framework  has  been  developed  in  the  literature.  A  methodology  needs 
to  be  developed  to  obtain  confidence  bounds  on  the  residual  stress.  Robust  design  is  an  important 
aspect  that  must  be  addressed,  because  small  difference  in  estimation  of  compressive  residual 
stresses  leads  to  an  order  of  magnitude  in  the  fatigue  life  estimates. 

3.6  Project  Scope  and  Organization 

This  section  describes  the  scope  of  the  project.  Ligure  10  provides  the  overview  of  the  project. 
The  main  goal  of  the  research  is  to  develop  a  method  to  quantify  uncertainty  in  the  LP  residual 
stresses  due  to  regression  uncertainty. 
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Figure  10:  Project  Scope 


The  framework  is  comprised  of  three  components:  (i)  regression  uncertainty  quantification  for 
material  model,  (ii)  uncertainty  propagation  to  LP  residual  stress,  and  (iii)  confidence  band  on  LP 
residual  stresses.  Regression  uncertainty  represents  the  variability  in  the  model  constant  estimates 
from  the  non-linear  regression  analysis.  A  statistical  technique  known  as  bootsrap  for  regression 
is  adapted  to  validate  the  normality  assumption  of  the  model  constant  estimates.  For  the 
uncertainty  propagation,  a  Taylor  series  expansion  is  used  as  an  approximation  to  propagate  the 
regression  uncertainty  through  the  FEA.  The  final  confidence  band  of  the  LP  residual  stress  is 
obtained  for  the  entire  depth. 

LP  has  many  process  parameters  that  must  be  precisely  determined  for  an  accurate  finite  element 
simulation.  Among  them  are  pressure  profile,  laser  spot,  and  material  model.  This  research 
investigates  different  material  models  and  validates  them  with  experimental  results.  An 
optimization-based  methodology  is  implemented  for  identifying  material  models  to  simulate  the 
LP  residual  stresses  when  very  little  experimental  data  of  material  behavior  is  available.  The 
project  is  organized  as  follows. 

Chapter  4  describes  the  details  of  the  finite  element  simulation  procedure.  Ligure  1 1  shows  the 
four  major  components  of  an  LP  simulation.  In  this  work,  research  is  conducted  in  the  material 
behavior  component.  Material  model  validation  for  the  simulation  of  the  LP  residual  stress  is 
described. 
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Figure  11:  Simulation  Overview 


Three  different  constitutive  material  models  are  used  for  validation  purposes:  the  EPP  model,  the 
JC  model  and  the  ZA  model.  The  analysis  procedure  and  geometric  modeling  are  also  described. 

In  Chapter  5,  simulation  of  the  LP  process  for  Inconel®718  material  is  described.  An 
optimization-based  methodology  is  implemented  to  predict  the  LP  residual  stresses.  This  inverse 
optimization  approach  obtains  the  model  constants  from  the  experimental  data  at  one  peak 
pressure  and  predicts  the  residual  stress  at  other  peak  pressures.  The  advantage  of  this  approach  is 
that  the  results  obtained  from  the  simulation  will  be  consistent  with  the  experimental  results.  The 
power  of  the  approach  can  be  seen  when  there  is  very  little  or  no  stress-strain  data  available  at 
different  strain  rates. 

In  Chapter  6,  the  framework  to  obtain  the  confidence  bounds  on  the  residual  stress  is  described.  A 
generic  definition  of  the  bootstrapping  technique  is  described  and  an  example  is  provided.  The 
technique  of  bootstrapping  for  regression  is  also  described.  The  two  categories  of  bootstrapping 
for  regression  are  defined,  and  an  example  for  each  method  is  provided.  Two  demonstration 
examples  are  provided  to  validate  the  uncertainty  framework.  The  developed  framework  is 
implemented  for  the  JC  model  and  the  confidence  band  for  the  LP  residual  stresses  is  obtained. 
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4.0  MATERIAL  MODEL  VALIDATION  LOR  LASER  PEENED  RESIDUAL  STRESSES 


This  chapter  provides  brief  details  of  the  Laser  Peening  (LP)  process  with  reference  to  the 
simulation  procedure.  Constitutive  material  model  advantages  and  disadvantages.  The  simulation 
procedure  is  described  followed  by  results  and  discussion. 

4.1  Laser  Peening  Description  With  Reference  to  Simulation 

The  LP  process  is  comprised  of  two  parts.  The  first  part  is  the  quick  process  of  laser  firing  and  the 
second  part  is  the  slow  recovery  of  the  component  to  equilibrium.  The  laser  firing  part  is  a  highly 
intense  process.  A  high  powered  laser  with  intensity  of  the  order  of  1 09  W /cm1  is  focused  on  an 
opaque  overlay.  The  laser  is  in  contact  with  the  opaque  overlay  for  a  few  nano  seconds.  Firing  the 
laser  on  the  component  creates  shock  waves  that  travel  inside  the  component  to  create 
compressive  residual  stresses  on  the  surface  regions.  The  relaxation  of  the  material  to  achieve 
equilibrium  is  a  relatively  slow  process  compared  to  laser  firing.  During  the  simulation  of  the  LP 
process,  this  mismatch  of  time  is  modeled  using  ABAQUS/Explicit  and  ABAQUS/Standard  for 
efficient  computation. 

4.1.1  Thick  and  Thin  Geometries 

This  section  describes  the  overview  of  the  simulation  models  that  can  be  developed  for  the  LP 
simulation.  Figure  12  shows  the  overview. 
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Simulation  Models 


Figure  12:  Overview  of  Simulation  Models 


The  simulation  models  can  be  divided  into  two  categories,  i)  thick  specimen  models  ii)  thin 
specimen  models.  Thick  specimen  models  are  defined  as  the  models  where  the  shock  waves  that 
are  reflected  from  the  bottom  surface  are  not  significant.  There  is  no  standard  way  of  defining  the 
exact  thickness  that  differentiates  the  specimen.  It  depends  on  the  material  being  used  and  also  on 
the  process  parameters.  Infinite  elements  are  used  for  the  thick  specimens  that  represent  the 
boundary  condition  of  shock  waves  not  reflecting  back.  Isotropic  hardening  works  well  for  the 
thick  specimens.  For  modeling  the  thin  specimens,  infinite  elements  cannot  be  used  since  the 
shock  waves  that  reflect  back  are  significant  in  prediction  of  the  residual  stresses.  A  combined 
hardening  model  must  be  used  since  isotropic  hardening  does  not  capture  the  effect  of  shock 
waves  from  both  the  direction.  All  the  models  developed  in  this  work  are  for  thick  geometries. 

4.2  Material  Model  Details 

The  plastic  deformation  during  the  LP  process  is  caused  by  laser- induced  shock  waves.  The 
waves  that  causes  stresses  beyond  the  yield  limit  is  generally  referred  to  as  plastic  waves.  These 
plastic  waves  are  known  as  shock  waves  when  the  wave  front  is  steep.  LP  is  a  very  high  strain  rate 
process.  The  component  can  experience  strain  rates  up  to  106/v,  which  is  very  high  compared  to 
conventional  strain  rates.  At  low  strain  rates  the  material  behavior  is  independent  of  strain  rates. 
The  material  response  to  impact  loading  may  be  significantly  different  from  static  or  quasi-static 
loading  [68].  Accurate  material  modeling  is  needed  to  simulate  the  material  behavior.  Various 
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constitutive  models  have  been  suggested  by  researchers  such  as  Bodner  [69],  Miller  [70], 
Bommann  et  al  [71],  Johnson  and  Cook  [72]  and  Zerilli  and  Armstrong  [73]for  different 
applications.  The  later  two  models  are  investigated  in  this  work  along  with  the  EPP  model 
because  of  their  popularity  in  literature. 

4.2.1  The  EPP  Model 


In  the  EPP  model,  no  strain  hardening  and/or  strain  rate  dependence  of  flow  stress  is  considered. 
Once  the  plastic  regime  is  reached,  the  stress  remains  constant.  The  yield  stress  is  derived  based 
on  Hugoniot  elastic  limit  (HEL)  because  LP  is  a  shock  wave  phenomenon.  HEL  is  defined  as  the 
axial  stress  required  for  plastic  deformation  under  uniaxial  strain  conditions.  It  is  assumed  that 
yielding  occurs  when  the  stress  in  the  direction  of  shock  wave  reaches  the  HEL.  The  relationship 
between  HEL  and  yield  strength  (oy)  is  shown  in  Equation  1: 


(1) 


where  v  is  Poisson’s  ratio.  The  advantage  of  this  model  is  that  only  a  small  amount  of  data  is 
required  to  estimate  the  yield  stress.  The  disadvantage  is  that  it  does  not  account  for  strain 
hardening  and/or  strain  rate  dependence. 


4.2.2  The  JC  Model 


The  JC  model  is  one  of  the  most  frequently  used  models  for  impact  studies  [72].  The  JC  model 
describes  the  flow  stress  of  the  material  as  a  product  of  three  terms:  a  strain  hardening  term,  a 
strain  rate  dependent  term,  and  a  temperature  term.  It  is  described  in  Equation  2: 


g  =  [A  +  Be”] 


£ 

/  " 

1  +C//7  — 

\-Tm 

Co 

- 

(2) 


where  A,  B ,  n,  C,  £q,  and  m  are  experimentally  determined  constants.  £  and  £  are  strain  and  strain 
rate,  respectively.  T  is  a  non-dimensionalized  temperature: 


T  = 


T-Tr 
/  —  Tm 


(3) 


Tr  and  Tm  are  room  temperature  and  melting  temperature,  respectively.  Parameter  A  is  the  initial 
yield  strength  at  room  temperature.  Parameters  B  and  n  take  strain  hardening  into  account,  while 
parameter  m  models  the  thermal  softening.  Parameter  C  represents  strain  rate  sensitivity.  The 
main  advantage  of  the  JC  model  is  that  the  estimation  of  the  parameters  is  simple  and  easy 
because  it  allows  for  isolation  of  the  three  effects. 
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4.2.3  The  ZA  Model 


The  ZA  constitutive  model  is  based  on  dislocation  mechanics  and  the  crystal  structure  of  the 
material  [73].  There  are  several  generations  of  the  model.  Initially,  the  model  addressed  Face 
Centered  Cubic  (FCC)  and  Body  Centered  Cubic  (BCC)  structures.  The  flow  stress  relationship  is 
shown  in  Equations  4  and  5  for  FCC  and  BCC  structures,  respectively: 

a  =  Ci+C5enexp  {-C3T +  C4Tlne)  (4) 

<j  =  Ci  +C2exp(— C3T  +  C4Tln£)  -\-C3£n  (5) 

Cj,  C2,  C3,  C4,  C5,  and  n  are  material  constants  that  need  to  be  determined.  Ci,  C5,  and  n  are 
similar  to  A,  B.  and  n,  respectively,  of  the  JC  model.  The  work  has  been  extended  to  include  strain 
recovery  for  Hexagonal  Closely  Packed  (HCP)  metals  [74].  Equation  6  shows  the  constitutive 
model 

a  =  oa  +  Be~^^Int)T  +B0  sj er  (l  -  e~el^)e~^~ailnt)T  (6) 

where  oa,  B.  /30,  f5\ ,  Bq,  0Cq,  and  oqare  the  constants  that  need  to  be  estimated.  The  term  in  the 
square  root  is  the  strain  recovery  term.  This  addresses  the  shear  instability,  an  important 
consideration  for  HCP  structures  (Ti-6A1-4V).  The  ZA  model  considers  the  interaction  effect 
between  strain  rate  and  temperature.  The  disadvantage  of  the  ZA  model  is  that  it  has  a  high 
number  of  constants  that  must  be  determined. 

4.2.4  Parameter  Estimation 

Ti-6A1-4V  material  is  used  in  this  work.  The  material  properties  are  listed  in  Table  1  [38]. 


Table  1:  Material  Properties  of  Ti-6A1-4V 


Property 

Value 

Young’s  Modulus  (GPa) 

113.8 

Poisson’s  Ratio 

0.342 

Density  (kgm  3) 

4500 

HEL  (MPa) 

2800 

The  material  constants  for  the  JC  and  ZA  models  are  estimated  by  fitting  the  experimental  data  of 
stress-strain  curves  for  different  strain  rates  [75].  Although  more  recent  experimental  data  is 
available  in  the  literature  [76],  these  results  are  not  used  because  the  experiments  were  conducted 
on  low  grade  Ti-6A1-4V,  which  is  not  used  for  aerospace  applications.  Higher  strain  rate 
experiment  data  [77]  was  found  and  will  be  used  for  the  simulation  in  the  next  chapter.  Equation 
6,  which  considers  the  strain  recovery  term,  is  used  as  the  ZA  material  model.  Temperature 
effects  are  not  considered  in  this  work  because  LP  is  a  mechanical  process  [78].  MATLAB’s 
lsqcurvefit  [79]  function  is  used  for  fitting  the  experiment  data,  lsqcurvefit  is  a  nonlinear  curve 
fitting  function  in  a  least  square  sense.  Figures  13  and  14  show  the  fit  to  the  experimental  data  for 
the  JC  and  ZA  models  respectively. 
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Figure  13:  The  JC  Model  Fit  to  Experimental  Data  [75]  of  Multiple  Strain  Rates 
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Flow  Stress[MPa] 


Figure  14:  The  ZA  Model  Fit  to  Experimental  Data  [75]  of  Multiple  Strain  Rates 


The  constants  obtained  for  the  JC  model  and  the  ZA  model  are  shown  in  tabular  format  in  Tables 
2  and  3. 


Table  2:  Model  Constants  for  the  JC  Model 


Parameter 

Value 

A 

950.228  MPa 

B 

603.3825  MPa 

n 

0.1992 

C 

0.0198 

£o 

0.0009/s 
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Flow  Stress  [MPa] 


Table  3:  Model  constants  for  the  ZA  Model 


Parameter 

Value 

<ya 

945. 1961  MPa 

B 

246.6467MPa 

Po 

1.1636  x  1CT6 

A 

0.1065 

Bo 

1481.249  MPa 

£r 

0.0538 

CC 0 

Hr5 

a\ 

3.1564  x  10^4 

Figures  15  and  16  show  the  strain  rate  dependence  curves  at  different  strains  for  the  JC  and  ZA 
models. 


Figure  15:  Strain  Rate  Dependence  Curve  for  the  JC  Model 
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Figure  16:  Strain  Rate  Dependence  Curve  for  the  ZA  Model 


These  figures  indicate  that  for  the  JC  model,  the  stress  varies  linearly  with  strain  rate  and  the 
curves  are  approximately  parallel  for  different  strains.  Nonlinear  variation  is  observed  for  the  ZA 
model.  In  the  ZA  model,  more  strain  rate  hardening  is  observed  for  higher  strain.  This  model 
could  be  closer  to  the  real  behavior  because  the  material  hardens  at  a  faster  rate  at  higher  strains. 

4.3  Simulation  Procedure 

In  this  work,  FEA  is  used  to  simulate  the  LP  process.  ABAQUS  [80]  software  is  used  to  perform 
FEA.  Figure  17  shows  the  two-dimensional  representative  mesh.  It  is  an  axi-symmetric  model 
made  up  of  finer  and  infinite  elements.  CAX4R  (continuum  axi-symmetric  4  noded  reduced 
integration)  elements  are  used  for  finer  elements  and  CINAX4  (continuum  infinite  axi-symmetric 
4  noded)  elements  are  used  for  the  infinite  part.  LP  is  a  highly  localized  process. 


25 

Approved  for  public  release;  distribution  unlimited. 


Pressure  Loading 


Figure  17:  Representative  Axi-symmetric  FE  Mesh 


The  size  of  the  region  treated  with  the  shot  is  very  small  relative  to  the  component  size.  Finer 
elements  are  modeled  for  the  treated  region  and  infinite  elements  represent  the  rest  of  the 
component.  The  dimension  of  the  finer  mesh  region  depends  on  the  size  of  the  laser  spot. 
Researchers  [37,  41]  have  used  2  times  and  3.5  times  the  laser  spot  size  for  the  finer  mesh.  In  this 
work,  the  radius  of  the  laser  spot  is  2.8mm,  and  the  dimensions  of  the  fine  mesh  region  are  6mm  x 
6mm,  with  symmetric  boundary  conditions  as  shown  in  Figure  17.  The  results  with  higher 
dimensions  for  finer  mesh  are  very  similar  to  the  current  mesh  of  6mm  x  6mm.  To  achieve 
computational  efficiency,  the  dimensions  of  6  mm  x  6  mm  were  used.  A  refined  mesh  is  required 
to  accurately  capture  the  shock  wave  imparted  by  the  laser.  A  mesh  size  of  241  x  241  nodes  is 
used  for  finer  elements  after  performing  the  convergence  study.  The  element  size  is  0.025  mm  x 
0.025  mm  for  the  mesh  size  used  for  the  finer  mesh. 

4.3.1  Pressure  Loads 

Pressure  loading  is  applied  on  the  top  surface  of  the  model  as  shown  in  Figure  17.  The  shock 
wave  generated  during  the  LP  process  is  applied  as  a  pressure  pulse  in  FEA.  There  are  different 
forms  of  pulse  shapes  used  in  the  literature.  The  shape  of  the  pressure  pulse  depends  on  the  laser 
system.  Ding  and  Ye  [78]  used  a  Gaussian  shape  and  Braisted  and  Brockman  [37]  used  a 
triangular  shape.  The  pressure  pulse  applied  in  this  work  is  shown  in  Figures  18  and  19  and  is 
similar  to  the  one  used  by  Nam  [38].  It  has  both  time  and  spatial  variation.  A  plot  of  peak 
pressure  vs.  power  density  was  developed  using  quartz  gages  and  this  plot  was  used  to  estimate 
the  average  peak  pressure  profile.  The  pressure  rises  for  the  first  few  nanoseconds  (ns)  and 
gradually  decreases  after  that. 
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Figure  18:  Temporal  Pressure  Profile  of  Shock  Wave 


Figure  19:  Spatial  Pressure  Profile  of  Shock  Wave 


4.3.2  Material  Models 

For  all  the  FEA  performed  in  this  work,  finer  elements  are  assumed  to  be  elastic-plastic  in  nature, 
while  the  infinite  elements  are  assumed  to  be  elastic  only  because  the  shock  wave  affects  a  small 
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portion  of  the  component.  These  distinctions  help  to  achieve  efficient  computation.  As  mentioned 
previously,  the  EPP,  JC  and  ZA  material  models  are  investigated. 


4.3.3  Analysis  Procedure 

The  analysis  can  be  categorized  into  two  stages:  dynamic  loading  analysis  and  static  equilibrium 
analysis.  During  the  first  stage,  the  high  speed,  intense  pressure  transient  loading  is  modeled  until 
all  plastic  deformation  has  taken  place.  The  second  stage  is  the  static  equilibrium  analysis  that  is 
performed  to  obtain  the  residual  stress  field.  The  data  from  the  end  of  the  dynamic  analysis,  such 
as  nodal  stresses,  strains,  and  displacements  are  transferred  to  the  equilibrium  analysis  using  a 
restart  file  in  ABAQUS.  After  the  static  equilibrium  analysis  is  completed,  there  are  two  options: 
(i)  one  more  LP  shot  or  (ii)  no  further  LP  shots.  If  no  further  LP  shots  are  required,  residual  stress 
data  is  obtained.  If  further  LP  shots  are  required,  the  residual  stress  data  obtained  from 
equilibrium  analysis  act  as  a  starting  point  for  the  new  LP  shot.  The  procedure  is  shown 
schematically  in  Ligure  20. 


28 

Approved  for  public  release;  distribution  unlimited. 


Figure  20:  Schematic  Representation  of  Analysis  Procedure 


4.3.3. 1  Dynamic  Loading  Analysis 

ABAQUS/Explicit  is  used  for  the  loading  part  of  the  analysis.  Two  of  the  important  parameters  to 
be  determined  for  the  analysis  are  (i)  time  increment  and  (ii)  final  time  to  complete  the  analysis. 
Time  increment  (At)  plays  an  important  role  in  the  convergence  and  accuracy  of  the  analysis.  If 
the  time  increment  is  greater  than  the  critical  limit  (A tcrit),  it  may  lead  to  an  unbounded  solution 
[80].  There  are  different  ways  to  estimate  the  critical  limit  [78].  A  simple  estimate  is  based  on 
element  length  and  wave  speed  of  the  material  (Cj).  The  critical  limit  is  calculated  as  Atcnf  =  ^r, 

where  L  is  the  length  of  the  smallest  element  and  Cj  —  \JE  /  p.  where  E  and  p  are 
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Young ’modulus  and  density,  respectively.  The  time  steps  for  each  model  are  optimized  from 
convergence  studies.  The  loading  analysis  is  performed  until  all  the  plastic  deformation  has  taken 
place.  This  is  shown  in  Figure  21,  the  history  plot  of  kinetic  energy.  This  plot  indicates  that  the 
final  time  of  4000  ns  is  suitable  for  LP. 


Figure  21:  History  of  Kinetic  Energy 


4.3.3.2  Static  Equilibrium  Analysis 

ABAQUS/Standard  is  used  for  the  equilibrium  analysis.  The  state  of  the  deformed  component 
from  the  end  of  dynamic  loading  analysis  is  transferred  using  a  restart  file  to  achieve  static 
equilibrium.  This  step  is  used  to  calculate  the  residual  stress  profile. 

4.3.3.3  Comparison  between  Explicit  and  Implicit  Algorithms 

ABAQUS/Explicit  and  ABAQUS/Standard  are  the  two  solvers  used  by  ABAQUS  in  the  FEA 
procedure.  ABAQUS/Explicit  uses  an  explicit  algorithm  and  ABAQUS/Standard  uses  implicit 
algorithm  for  the  numerical  analysis.  This  section  describes  a  brief  comparison  between  the  two 
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algorithms.  Explicit  and  implicit  are  direct  integration  methods  for  solving  the  equation  of  motion 
(Equation  7) 

[M]  {D}  +  [C]  {D}  +  [K]  {£>}  =  {Rext}  (7) 

where  [M\,  [C],  and  [K]  are  mass,  damping,  and  stiffness  matrices  respectively.  { D }  is  the 
displacement  vector.  { Rext }  corresponds  to  external  loads.  These  loads  are  in  general  a  function 
of  time.  The  direct  integration  methods  are  alternative  to  modal  analysis  methods  [81].  In  direct 
integration  methods,  displacement  vector  is  associated  with  subscript  n  to  denote  the  time  step. 
The  solution  methodology  for  explicit  algorithm  have  the  form  shown  in  Equation  8: 

{»}„+!  =  /  ({£>}„ ,  M„ ,  {D}„ ,  {£>}„_! , ...)  (8) 

This  displacement  {D}n+1  is  determined  in  terms  of  completely  historical  information  of 
displacements  and  its  time  derivatives.  The  displacement  at  the  current  step  is  a  function  /  of 
displacements  and  their  time  derivatives  of  previous  steps.  Implicit  methods  have  the  form  shown 
in  Equation  9: 

{»}„+!  (9) 

The  computation  of  {D}n+l  requires  knowledge  of  the  time  derivatives  of  {D}n+l.  The  implicit 
algorithm  determines  the  solution  with  iterations,  but  the  explicit  algorithm  determines  the 
solution  by  advancing  the  kinematic  state  from  the  previous  state.  Central  difference  scheme  is 
popularly  used  for  the  explicit  algorithm.  In  comparison,  the  explicit  algorithm  is  a  clear  choice 
for  wave  propagation  analysis,  especially  for  short  duration.  The  implicit  algorithm  is  better 
suited  for  long  duration  to  achieve  equilibrium. 

4.4  Results  and  Discussion 

Three  peak  pressures  of  5.5,  6.1,  and  8.3GPa  are  considered  in  this  work.  These  three  peak 
pressures  are  selected  based  on  the  availability  of  experimental  data  [38].  Three  LP  shots  are  used 
for  a  peak  pressure  of  6.1GPa,  and  one  LP  shot  is  used  for  the  other  two  peak  pressures.  The 
in-depth  experimental  results  of  radial  residual  stresses  at  a  distance  of  1mm  from  the  center  for 
all  three  peak  pressures  are  shown  in  Figures  22,  23,  and  24  for  the  EPP  model,  the  JC  model,  and 
the  ZA  model  respectively. 
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Figure  22:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  EPP 
model 
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Figure  23:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  JC 
model 
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Depth  [mm]  at  distance  of  1mm  from  center 


Figure  24:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  ZA 
model 

The  x-ray  diffraction  technique  was  used  to  measure  the  residual  stresses  [38].  This  is  a 
destructive  technique  used  to  measure  residual  stresses  along  depth.  First,  the  measurement  is 
made  on  the  surface.  A  thin  layer  is  removed  using  electropolishing  and  the  procedure  is  repeated 
for  each  depth  measurement.  This  technique  averages  the  residual  stress  over  a  square  area.  To 
compare  the  results  of  FEA  and  experiments,  the  same  kind  of  averaging  is  performed  for  stresses 
resulting  from  FEA.  A  1mm  x  1mm  area  is  considered,  and  the  x  and  y  coordinates  of  the 
locations  are  converted  into  cylindrical  coordinates  as  r  —  \Jx2  +  y2and  0  =  Tan  1  [y /x\ .  The 
transformation  of  the  stress  components  from  cylindrical  to  rectangular  is 
ox  —  orcos29  +  Oosin29  —  2z,Qsin9cos9 .  The  shear  component  (zrg)  is  zero  (definition  of 
axi- symmetric).  The  results  comparison  is  discussed  for  each  model  separately. 

4.4.1  EPP  Model 

Figure  25  shows  the  equivalent  plastic  strain  variation  versus  depth  for  different  peak  pressures. 
The  plastic  strain  data  is  obtained  at  a  distance  of  1mm  from  the  center  of  the  spot.  The  strains  are 
higher  on  the  surface  and  decrease  with  depth.  The  higher  the  magnitude  of  peak  pressure,  the 
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greater  is  the  magnitude  of  plastic  strain.  This  shows  that  the  model  is  predicting  the  general 
trends  with  respect  to  strain. 


Depth[mm]  at  a  distance  of  1  mm  from  center 

Figure  25:  Equivalent  Plastic  Strain  Variation  with  Depth  for  the  EPP  Model 


Figure  22  shows  the  radial  residual  stress  variation  with  depth,  for  different  peak  pressures.  It 
includes  both  simulation  and  experimental  results.  A  comparison  between  the  simulation  results 
for  different  peak  pressures  indicates  that  a  higher  magnitude  of  compressive  residual  stresses  are 
generated  for  greater  peak  pressures  only  after  a  certain  depth.  For  surface  and  near-surface 
regions,  the  model  is  not  able  to  predict  the  experimentally  observed  trends.  A  higher  magnitude 
of  stress  is  obtained  for  lower  peak  pressure,  and  for  all  three  peak  pressures,  there  is  a  decrease 
in  the  magnitude  of  residual  stress  on  the  surface  and  near-surface  region.  The  propagation  of 
stress  waves  into  the  material  for  peak  pressures  of  5.5  GPa  for  the  EPP  model  is  shown  in  Figure 
26  and  for  the  8.3GPa  is  shown  in  Figure  27.  Figures  28  and  29  show  the  stress  wave  propagation 
of  peak  pressures  of  5.5  GPa  and  8.3  GPa  for  the  JC  model  respectively. 
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Figure  26:  Propagation  of  Stress  Waves  for  the  EPP  model  at  5.5  GPa  Peak  Pressure 
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Figure  27:  Propagation  of  Stress  Waves  for  the  EPP  model  at  8.3  GPa  Peak  Pressure 
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Figure  28:  Propagation  of  Stress  Waves  for  the  JC  Model  at  5.5  GPa  Peak  Pressure 
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Figure  29:  Propagation  of  Stress  Waves  for  the  JC  model  at  8.3  GPa  Peak  Pressure 


During  the  first  200  ns  the  compressive  wave  declines  steeply  at  the  surface  for  both  peak 
pressures  for  the  EPP  model.  This  decline  is  predominant  for  a  peak  pressure  of  8.3GPa.  There  is 
no  decline  observed  in  the  stress  wave  propagation  for  the  JC  model  for  either  of  the  peak 
pressures.  This  initial  stress  wave  has  a  significant  effect  on  the  residual  stress,  due  to  the 
magnitude  of  the  wave.  A  comparison  between  the  simulation  and  the  experimental  data  for  three 
peak  pressures  of  5.5,  6.1  and  8.3GPa  can  be  seen  in  Figure  30. 
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Figure  30:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  EPP  Model 


For  a  peak  pressure  of  5.5GPa,  the  EPP  model  overestimates  the  residual  stresses  relative  to  the 
experimental  data  for  the  entire  depth  of  2mm.  The  EPP  model  overestimates  the  stresses  after  a 
certain  depth  for  peak  pressures  of  6.1  and  8.3GPa.  The  EPP  model  is  generally  considered  to  be 
a  conservative  model  for  prediction  of  internal  stresses.  However,  for  residual  stress 
determination,  the  EPP  model  always  overestimates  compared  to  experimental  data  since  the  EPP 
model  does  not  take  strain  hardening  into  consideration.  In  the  LP  process,  strain  rate  plays  an 
important  role,  and  the  current  EPP  model  does  not  consider  strain  rate  dependence. 

4.4.2  JC  Model 

Figure  31  shows  the  equivalent  plastic  strain  variation  with  depth,  for  different  peak  pressures. 
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Figure  31:  Equivalent  Plastic  Strain  Variation  with  Depth  for  the  JC  Model 


The  same  trends  as  observed  in  the  EPP  model  can  be  seen  here.  There  is  a  higher  magnitude  of 
plastic  strain  obtained  for  greater  peak  pressure.  However,  the  magnitude  of  plastic  strain  is  much 
lower  than  the  EPP  model  for  the  three  peak  pressures.  For  the  peak  pressure  of  8.3  GPa,  a  steep 
decrease  in  the  strain  can  be  seen  up  to  a  depth  of  0.5mm,  then  a  gradual  decrease  is  observed.  A 
comparison  between  radial  residual  stresses  for  three  peak  pressures  is  shown  in  Figure  32. 
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Depth  [mm]  at  distance  of  1  mm  from  center 


Figure  32:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  JC  Model 

Figure  32  indicates  that  the  JC  model  is  able  to  predict  the  trends  for  the  entire  depth  of  2  mm, 
unlike  the  EPP  model.  A  higher  magnitude  of  stress  is  obtained  for  greater  peak  pressure.  The 
magnitudes  of  stress  at  the  surface  for  peak  pressures  of  6.1  and  8.3  GPa  are  nearly  equal.  This 
can  be  attributed  to  the  number  of  shots.  For  the  peak  pressure  of  6.1GPa,  the  number  of  shots  is 
three,  while  for  the  peak  pressure  of  8.3  GPa,  it  is  only  one.  A  comparison  between  the  simulation 
and  the  experimental  data  at  all  three  peak  pressures  for  the  JC  model  can  be  seen  in  Figure  32. 
The  JC  model  is  in  better  agreement  with  experimental  results  for  the  three  peak  pressures, 
indicating  the  consistency  of  the  model.  However,  the  JC  model  overestimates  the  residual  stress 
on  the  surface  compared  with  the  experimental  data  for  peak  pressures  of  5.5  and  6.1  GPa.  The 
JC  model  maintains  the  trend  of  overestimating  the  residual  stress  compared  with  experimental 
data  for  the  three  peak  pressures  (except  from  two  data  points  for  peak  pressure  of  6.  lGPa).  In  the 
next  section,  the  results  for  the  ZA  model  are  discussed.  Unlike  the  JC  model,  the  ZA  model 
overestimates  on  the  surface  and  underestimates  at  higher  depth. 
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4.4.3  ZA  Model 


Figure  33  shows  the  equivalent  plastic  strain  variation  with  depth  for  the  three  peak  pressures. 


Depthfmm]  at  a  distance  of  1mm  from  center 


Figure  33:  Equivalent  Plastic  Strain  Variation  with  Depth  for  the  ZA  Model 


The  ZA  model  is  able  to  predict  the  general  trend  of  higher  magnitude  of  plastic  strain  for  greater 
peak  pressure.  The  magnitude  of  strain  predicted  by  the  ZA  model  at  the  near-surface  region  is  in 
between  the  predictions  of  the  EPP  and  JC  models.  However,  at  higher  depth,  the  ZA  model  has 
the  lowest  strain,  followed  by  the  EPP  model,  and  the  JC  model  as  the  highest  magnitude  of 
strain.  There  is  also  a  steep  decrease  in  strain  at  the  near-surface  region  as  seen  in  the  other  two 
models.  The  plastic  strain  curves  at  peak  pressures  of  5.5  and  6.1GPa  are  relatively  close  to  each 
other  compared  with  the  other  two  models.  Radial  residual  stress  variation  with  depth  for 
different  peak  pressures  is  shown  in  Figure  34. 
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Depth  [mm]  at  distance  of  1mm  from  center 


Figure  34:  Residual  Stress  Comparison  Between  Simulation  and  Experiment  for  the  ZA  Model 


The  ZA  model  is  also  able  to  predict  the  trend  for  the  entire  depth.  A  higher  magnitude  of  stress  is 
obtained  for  greater  peak  pressures.  The  magnitude  of  stress  at  the  surface  for  peak  pressures  of 
6.1  and  8.3  GPa  are  nearly  equal,  as  seen  in  the  JC  model.  The  magnitude  of  stress  at  the  surface 
and  near-surface  regions  is  high  compared  to  the  JC  model.  A  comparison  between  the  simulation 
results  and  the  experimental  data  shows  that  the  ZA  model  overestimates  the  residual  stress  on  the 
surface  and  near-surface  regions  for  all  three  peak  pressures.  On  the  surface  it  overestimates  not 
only  with  respect  to  the  experimental  data,  but  also  with  respect  to  the  JC  and  EPP  models.  At 
higher  depths,  even  though  the  ZA  model  underestimates  relative  to  the  experimental  data,  there 
seems  to  be  a  better  agreement  with  experimental  results  of  residual  stress  results  on  the  surface. 

4.4.4  Additional  performance  metrics 

For  further  comparison  between  the  models,  four  performance  metrics  are  developed:  (a) 
maximum  compressive  stress,  (b)  maximum  tensile  stress,  (c)  area  of  compressive  stress  and  (d) 
depth  of  compressive  stress.  The  area  of  compressive  stress  is  a  normalized  area.  It  is  calculated 
as  the  ratio  of  the  number  of  nodes  with  compressive  stress  to  the  total  number  of  nodes.  The 
depth  of  the  compressive  stress  is  considered  at  a  distance  of  1mm  from  the  center.  A  region  with 
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a  depth  of  2  mm  is  considered,  and  a  distance  of  0.5mm  from  both  ends  is  left  out  for  the  above 
mentioned  metrics.  Tables  4,  5,  and  6  show  the  performance  metrics  for  peak  pressures  of  5.5 
GPa,  6.1  GPa,  and  8.3  GPa  respectively. 


Table  4:  Performance  Metrics  Comparison  Between  the  Models  for  Peak  Pressure  of  5.5  GPa 


Performance  Metrics 

Peak  Pressure  of  5.5  GPa 

EPP 

JC 

ZA 

Depth  of  Compressive  Stress  [mm] 

0.9959 

1.4938 

0.6473 

Area  of  Compressive  Stress  [normalized] 

0.5891 

0.8624 

0.5778 

Maximum  Compressive  Stress  [MPa] 

-669.84 

-600.60 

-778.38 

Maximum  Tensile  Stress  [MPa] 

90.24 

24.53 

27.8 

Table  5:  Performance  Metrics  Comparison  Between  the  Models  for  Peak  Pressure  of  6.1  GPa 


Performance  Metrics 

Peak  Pressure  of  6.1  GPa 

EPP 

JC 

ZA 

Depth  of  Compressive  Stress  [mm] 

1.2448 

1.668 

0.7967 

Area  of  Compressive  Stress  [normalized] 

0.6528 

0.9222 

0.6113 

Maximum  Compressive  Stress  [MPa] 

-745.84 

-738.48 

-968.62 

Maximum  Tensile  Stress  [MPa] 

112.40 

23.23 

37.53 

Table  6:  Performance  Metrics  Comparison  Between  the  Models  for  Peak  Pressure  of  8.3  GPa 


Performance  Metrics 

Peak  Pressure  of  8.3  GPa 

EPP 

JC 

ZA 

Depth  of  Compressive  Stress  [mm] 

1.8174 

1.9917 

1.4689 

Area  of  Compressive  Stress  [normalized] 

0.6528 

0.9925 

0.763 

Maximum  Compressive  Stress  [MPa] 

-642.20 

-778.12 

-963.21 

Maximum  Tensile  Stress  [MPa] 

57.05 

39.49 

63.51 

Compared  to  both  the  EPP  and  ZA  models,  the  JC  model  predicts  a  higher  depth  of  compression 
for  all  three  peak  pressures.  The  ZA  model  predicts  the  highest  compressive  stress  for  all  three 
peak  pressures,  while  the  JC  model  predicts  the  lowest  compressive  stress  for  peak  pressures  of 
5.5  and  6.1GPa.  The  JC  model  predicts  the  highest  area  of  compressive  stress  for  all  three  peak 
pressures,  which  shows  consistency  in  the  results. 
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4.5  Section  Summary 


Numerical  simulation  of  the  LP  process  is  addressed  in  this  chapter.  FEA  is  shown  to  be  a  useful 
method  to  simulate  the  process.  The  strain  rates  involved  in  the  LP  process  are  on  the  order  of 
10V\  and  experimental  results  of  material  behavior  are  typically  not  available  at  such  high 
strain  rates.  Accurate  material  models  help  in  providing  a  benchmark  simulation.  The  nonlinear 
regression  method  is  used  to  fit  the  experimental  data  to  obtain  the  parameters  for  the 
corresponding  models.  The  EPP  model,  which  is  most  often  used  in  the  literature,  is  shown  to 
produce  inconsistent  results.  The  ZA  model,  which  is  based  on  dislocation  mechanics,  produces 
consistent  trends  but  overestimates  the  results  compared  to  experimental  data.  The  JC  model  is 
shown  to  produce  consistent  results  matching  the  trends  and  better  agreement  with  experimental 
results. 
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5.0  INVERSE  OPTIMIZATION  OF  MATERIAL  MODELS  FOR  SIMULATION  OF 

LASER  PEENED  RESIDUAL  STRESSES 


The  limitation  of  the  existing  simulation  methodology  is  described  in  this  chapter.  This  leads  to 
the  development  of  a  new  optimization-based  approach  for  simulation  of  residual  stresses  induced 
by  LP.  Brief  description  of  the  experimental  residual  stress  results  is  provided. 

5.1  Limitation  of  Existing  Simulation  Methodology 

The  existing  simulation  methodology  of  residual  stresses  induced  from  LP  is  based  on 
experimental  data  about  the  material  behavior.  Figure  35  shows  a  flow  chart  of  the  simulation 
process.  The  JC  model  is  used  as  an  illustration. 


Figure  35:  Traditional  Approach  for  Simulation  of  Laser  Peened  Residual  Stresses 


Stress-Strain  data  at  lower  strain  rates  of  a  material  is  used  to  curve  fit  the  material  model  under 
investigation.  The  material  behavior  is  extrapolated  at  higher  strain  rates  where  no  experimental 
data  is  available.  Stress-strain  data  at  different  strain  rates  is  typically  not  available  for  many 
materials.  Inconel®718  is  one  such  material.  The  limitation  of  the  existing  simulation 
methodology  demands  a  new  approach  to  simulate  the  residual  stresses  induced  by  LP. 
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5.2  Optimization-based  Approach 

The  proposed  optimization-based  approach  is  shown  in  Figure  36. 


Figure  36:  Optimization-Based  Approach  for  Simulation  of  Laser  Peened  Residual  Stresses 


Material  model  constants  are  the  design  variables  during  the  optimization  procedure.  An  initial 
guess  is  provided  for  the  material  model  constants.  These  initial  values  are  obtained  from  the  user 
experience  or  knowledge  about  the  material.  In  this  work,  the  initial  guess  for  Inconel®718  is 
based  on  user  experience  and  the  initial  guess  for  Ti-6A1-4V  material  is  based  on  results  from 
traditional  approach.  The  drawback  of  the  optimization-based  approach  is  that  the  optimum 
model  constants  obtained  are  dependent  on  the  initial  guess.  LP  simulation  is  performed  and  the 
compressive  residual  stress  is  computed  at  locations  for  which  the  experimental  data  is  available. 
The  simulated  compressive  residual  stress  is  compared  with  the  experimental  residual  stress  and 
the  error  between  them  is  quantified.  This  error  is  considered  as  the  objective  function.  An 
unconstrained  optimization  formulation  is  solved  to  obtain  the  optimum  model  constants  as 
shown  in  Equation  10: 
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Min  Error 

Lower  Bounds^  Model  Constants  <  Upper  Bounds 


(10) 

(ID 


The  error  is  minimized  and  the  model  constants  are  iteratively  solved  until  convergence.  The 
convergence  criterion  to  stop  the  iterations  is  that  the  difference  between  the  objective  function 
from  the  previous  iteration  should  be  less  than  a  threshold  value.  The  objective  function  used  in 
this  work  is  weighted  least  squares  as  shown  in  Equation  12: 


The  experimental  data  is  divided  into  three  regions  based  on  the  location.  w\,  w 2,  and  w 3  in 
Equation  12  act  as  weights  for  regions  I,  IT  and  III  respectively.  Details  of  each  region  are 
provided  in  the  results  section.  The  terms  RSexp  and  RSsjm  are  experimental  residual  stress  and  the 
residual  stress  predicted  by  the  simulation  respectively.  In  this  work  the  regions  are  weighted 
equally,  but  the  weights  can  be  changed  based  on  the  agreement  between  simulation  and 
experimental  results.  MATLAB’s  optimization  toolbox  [79],  that  uses  a  gradient-based  approach, 
is  used  to  perform  optimization.  Side  bounds  are  placed  considering  the  practicality  of  the  model 
constants.  The  above-mentioned  approach  is  performed  for  one  set  of  LP  conditions  and  optimum 
model  constants  are  obtained.  The  same  model  constants  are  used  for  other  operating  conditions 
and  the  residual  stress  predictions  are  compared  with  the  experimental  data  to  demonstrate  the 
consistency  of  the  approach.  The  advantage  of  this  approach  is  that  the  results  obtained  from  the 
simulation  will  be  consistent  with  the  experimental  results.  The  power  of  the  approach  can  be 
seen  when  there  is  very  little  or  no  stress-strain  data  available  at  different  strain  rates.  Using  this 
inverse  optimization  approach,  a  material  model  can  be  employed  to  simulate  the  residual  stresses 
induced  by  the  LP  process. 

5.3  Material  Model  Details 

Three  material  models  are  investigated  for  the  validation  of  the  optimization-based  method.  Two 
of  the  models  are  the  JC  and  the  ZA  models  that  were  previously  used.  Description  of  these  two 
models  are  provided  in  the  previous  chapter.  The  Khan-Huang-Liang  (KHL)  [77]  model  is  the 
third  material  model  under  investigation.  The  JC  and  KHL  models  are  considered 
phenomenological  models,  while  the  ZA  model  is  a  physics-based  model. 

5.3.1  KHL  Model 


The  KHL  model  is  based  on  the  JC  model.  The  flow  stress  equation  is  shown  in  Equation  13 


a 


A  +  B 


(13) 


A  major  feature  of  this  model,  compared  to  the  JC  model,  is  that  it  can  accommodate  decreasing 
work-hardening  with  an  increasing  strain  rate  through  the  constant  n\ .  £  and  e  are  the  plastic 
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strain  and  strain  rate,  respectively.  Model  constants  A  and  B  are  similar  to  the  constants  in  the  JC 
model.  A  represents  initial  yield  strength  at  room  temperature.  B  and  n  represent  strain  sensitive 
parameters.  Dq  is  an  arbitrarily  chosen  upper  bound  on  the  strain  rate.  In  this  work,  Dq  is 
assumed  to  be  108.  £*  is  the  reference  strain  rate  at  which  the  model  constants  are  determined.  In 
most  cases,  £*is  assumed  to  be  equal  to  lv_1 . 

5.4  Simulation  Procedure 

The  simulation  procedure  described  in  the  previous  chapter  is  used.  The  addition  of  KHL  model 
is  implemented  and  the  material  properties  [82]  for  Inconel®718  are  shown  in  Table  7. 


Table  7:  Material  Properties  for  Inconel®718 


Property 

Value 

Young’s  Modulus  (GPa) 

199.95 

Poisson’s  Ratio 

0.29 

Density  (kgm5) 

8193.25 

5.5  Results  and  Discussion 

For  Inconel®718,  the  residual  stress  predictions  through  depth  are  compared  with  experimental 
data  at  the  center  of  the  laser  spot,  while  the  comparison  is  performed  at  a  distance  of  1mm  from 
the  center  of  the  spot  for  Ti-6A1-4V.  The  averaging  procedure  described  in  the  previous  chapter  is 
implemented.  For  Ti-6A1-4V  material,  residual  stress  results  using  the  traditional  approach  are 
available  [83].  A  comparison  is  made  between  two  approaches  to  demonstrate  the  consistency  of 
the  optimization-based  approach.  The  residual  stress  results  for  Inconel®718  material  are 
presented  first,  followed  by  Ti-6A1-4V  material. 

5.5.1  Validation  of  Optimization-Based  Approach  for  Inconel®718 

Inconel®718  is  a  nickel-based  super  alloy  used  in  several  gas  turbine  components  that  operate  at 
high  temperatures  such  as  1 00(7' C.  The  experimental  procedure  to  perform  LP  experiments  is 
discussed  followed  by  comparison  with  simulation  predictions. 

5.5.1. 1  LP  Experimental  Procedure 

LP  experiments  were  performed  at  LSP  Technologies,  Dublin,  OH.  Four  coupons  were  laser 
peened  to  provide  experimental  residual  stress  profiles  to  compare  with  the  simulation  results. 
The  coupons  were  25  mm  x  25  mm  square  with  a  thickness  of  12.7  mm.  They  were  laser  peened 
with  one  laser  spot  with  a  diameter  of  5.6  mm  in  the  center  of  one  face.  The  surface  to  be  laser 
peened  was  electropolished  to  remove  machining  stresses  before  LP.  The  LP  was  performed  with 
an  Nd-glass  laser  using  black  paint  as  the  opaque  overlay  and  water  as  the  transparent  overlay. 
The  LP  conditions  are  shown  in  Table  8. 
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Table  8:  LP  Conditions  Used  for  Inconel®718 


Energy  Density  J/cm 2] 

Laser  Pulse  Length  [ns] 

Power  Density  [GW /cm2] 

105 

25 

4.2 

132 

21 

6.3 

181 

28 

7.2 

146 

18 

8.1 

X-ray  diffraction  method  was  used  for  the  measurement  of  residual  stress  profiles  from  the  laser 
peened  surface.  With  this  technique,  the  measurements  were  made  using  a  2  mm  square  x-ray 
diffraction  spot  on  the  coupon  surface,  centered  within  the  laser  spot.  After  taking  a  surface 
measurement,  the  material  was  incrementally  removed  from  surface  by  electropolishing  and 
another  measurement  was  taken  on  the  exposed  surface.  This  was  done  for  each  depth  in  the 
residual  stress  profile.  The  x-ray  diffraction  result  provides  a  measurement  of  the  elastic  strain. 
The  elastic  modulus  is  then  used  to  derive  the  residual  stress  at  each  depth.  Corrections  are 
applied  to  account  for  the  x-ray  penetration  depth  and  the  stress  relaxation  accompanying  the 
removal  of  material  from  the  surface.  The  residual  stress  data  at  different  depths  and  peak 
pressures  shown  in  Table  9  are  measured  at  LSP  Technologies  Inc,  provided  Dr.  Alan  H.  Clauer. 


5.5.1.2  Residual  Stress  Comparison  Between  Simulations  and  Experiments 

The  regions  are  divided  based  on  the  location  of  the  experimental  data.  Region  I  is  considered  up 
to  a  depth  of  0.5  mm  from  the  surface,  Region  II  is  from  0.5  mm  to  1.5  mm  and  Region  III  is  from 
1.5  mm  onwards.  The  optimization-based  approach  is  adapted  to  the  JC  model  and  the  KHL 
model.  A  peak  pressure  of  6.6  GPa  is  arbitrarily  chosen  to  perform  the  optimization.  The 
constants  obtained  for  the  JC  model  are  shown  in  Table  lO.The  model  constants  obtained  for  the 
KHL  model  are  shown  in  Table  16.  The  same  model  constants  are  used  for  other  peak  pressures 
of  5.3,  7.1,  and  7.6  GPa  to  demonstrate  the  consistency  of  the  inverse  optimization  approach. 
Figures  37,  38,  39,  and  40  show  the  residual  stress  comparison  between  simulations  and 
experimental  results  for  peak  pressures  of  5.3,  6.6,  7.1,  and  7.6  GPa  respectively. 
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Residual  Stress  [MPa] 


Table  10:  Model  Constants  Using  Optimization-Based  Approach  for  the  JC  Model  for  Inconel®718 


Parameter 

Value 

A 

506.44 MPa 

B 

1578.86MPi3 

n 

7.94  x  10~4 

C 

2.03  x  10~3 

£o 

6.45  x  10~4/s 

Figure  37:  Residual  Stress  Comparison  for  Inconel®718  with  5.3  GPa  Peak  Pressure 
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Depth  [mm] 


Figure  38:  Residual  Stress  Comparison  for  Inconel®718  with  6.6  GPa  Peak  Pressure 
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Figure  39:  Residual  Stress  Comparison  for  Inconel®718  with  7.1  GPa  Peak  Pressure 
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Depth  [mm] 


Figure  40:  Residual  Stress  Comparison  for  Inconel®718  with  7.6  GPa  Peak  Pressure 


Simulation  results  follow  the  trends  of  the  experimental  results  for  all  peak  pressures  in  both  the 
JC  and  the  KHL  models.  The  results  predicted  by  the  JC  model  and  the  KHL  model  are  similar  to 
each  other.  The  reason  for  this  could  be  attributed  to  the  fact  that  the  KHL  model  is  based  on  the 
JC  model.  The  simulation  results  show  a  consistent  trend  of  underestimating  the  experimental 
results  for  all  peak  pressures.  For  peak  pressures  of  5.3  and  6.6  GPa,  the  simulation  prediction 
performs  consistently  for  the  entire  depth.  For  peak  pressures  of  7.1  and  7.6  GPa,  the  simulation 
results  are  in  relatively  poor  agreement  with  the  experimental  data  as  compared  to  the  previous 
two  peak  pressures.  The  simulation  predictions  tend  to  overestimate  at  near  surface  regions  and 
underestimate  at  lower  depths.  The  error  between  the  experimental  data  and  simulation  results  is 
quantified.  Table  1 1  shows  the  least  square  error  between  the  experimental  data  and  simulation 
results  for  different  models.  The  error  is  computed  using  all  nine  experimental  data  locations  and 
the  corresponding  simulation  predictions.  Equation  10  is  used  to  obtain  the  least  square  error.  In 
Table  11,  the  lowest  error  is  for  peak  pressure  of  6.6  GPa  for  both  the  JC  and  the  KHL  models. 
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Table  11:  Least  Square  Error  Comparison  for  Inconel®718 


Inconel®7 1 8 

JC  Model 

KHL  Model 

Peak  Pressure  (GPa) 

5.3 

6.6 

7.1 

7.6 

5.5 

6.6 

7.1 

7.6 

Optimized  (MPa) 

261.76 

241.19 

547.94 

634.95 

288.82 

241.19 

553.04 

638.66 

Normalized 

1.08 

1.0 

2.27 

2.63 

1.19 

1.0 

2.29 

2.64 

The  error  is  normalized  to  the  least  error  value  (set  =  1 .0)  since  the  relative  error  is  a  more 
important  factor  than  the  absolute  value.  The  normalized  error  values  from  Table  1 1  demonstrates 
that  both  models  perform  very  similarly  at  all  peak  pressures.  Table  1 1  indicates  the  error 
between  the  experimental  data  and  simulation  result  at  a  peak  pressure  of  7.6GPa  is  2.64  times  the 
error  at  the  peak  pressures  of  5.3  GPa  for  the  KHL  model. 

5.5.2  Validation  of  Optimization-based  Approach  and  Comparison  with  Traditional 
Approach  for  Ti-6A1-4V 

Three  different  peak  pressures  of  5.5,  6.6,  and  8.3  GPa  were  selected  for  the  simulation  of  the 
residual  stresses  based  on  the  availability  of  the  experimental  data  [38].  Region  I  is  considered  up 
to  a  depth  of  0.2  mm  from  the  surface,  Region  II  is  from  0.2  mm  to  0.6  mm  and  Region  III  is  from 
0.6  mm  onwards.  The  optimization-based  approach  was  implemented  for  the  peak  pressure  of  5.5 
GPa  for  all  three  material  models  (the  JC,  the  ZA,  and  the  KHL  models).  The  same  model 
constants  are  used  for  the  other  two  peak  pressures  to  evaluate  the  consistency  of  the  approach. 
The  constants  obtained  for  the  JC  model  using  the  optimization-based  approach  are  shown  in 
Table  12: 


Table  12:  Material  Model  Constants  Using  Optimization-Based  Approach  for  the  JC  Model 


Parameter 

Value 

A 

1092.68Mftz 

B 

1266.94Mftz 

n 

0.0851 

C 

9.39  x  10  3 

£o 

9.0  x  10~4/s 

By  comparison,  the  model  constants  obtained  from  the  traditional  approach  [83]  are  shown  in 
Table  13 
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Table  13:  Material  Model  Constants  Using  Traditional  Approach  for  the  JC  Model  [83] 


Parameter 

Value 

A 

950.23 MPa 

B 

603.28 MPa 

n 

0.19921 

C 

0.0198 

£o 

9.32  x  10~7s 

The  constants  obtained  for  the  ZA  model  using  the  optimization-based  approach  are  shown  in 
Table  14 

Table  14:  Material  Model  Constants  Using  Optimization-Based  Approach  for  the  ZA  Model 


Parameter 

Value 

<Ja 

472. 89 A/Pa 

B 

247. 13 A/Pa 

/3o 

1.16  x  10  6 

Pi 

0.01598 

Bo 

2221. 87 MPa 

£r 

0.0308 

CCo 

Hr5 

ai 

1.85  x  10  4 

By  comparison,  the  model  constants  obtained  from  the  traditional  approach  are: 


Table  15:  Material  Model  Constants  Using  Traditional  Approach  for  the  ZA  Model  [83] 


Parameter 

Value 

<ya 

945.19 MPa 

B 

246. 65 MPa 

A) 

1.16  x  10  6 

Pi 

0.1065 

Bo 

1481.25 MPa 

£r 

0.0538 

CCo 

HP5 

3.16  x  10~4 

The  constants  obtained  for  the  KHL  model  using  the  optimization-based  approach  are  shown  in 
Table  16 
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Table  16:  Model  Constants  Using  Optimization-Based  Approach  for  the  KHL  Model 


Parameter 

Value 

A 

677.7 \MPa 

B 

9\1 .16\MPa 

n\ 

0.1526 

770 

0.0931 

c 

0.0415 

The  constants  obtained  for  the  KHL  model  using  the  traditional  approach  are: 


Table  17:  Material  Model  Constants  Using  Traditional  Approach  for  the  KHL  Model 


Parameter 

Value 

A 

1170.10M.Pa 

B 

812.511  MPa 

n\ 

0.4481 

no 

0.3301 

C 

0.0212 

Figure  41,  42,  and  43  shows  the  residual  stress  results  at  peak  pressures  of  5.5,  6.1,  and  8.3  GPa 
respectively  for  the  JC  model.  For  the  peak  pressure  of  5.5  GPa,  the  results  predicted  by  the 
optimization-based  approach  are  closer  to  the  experimental  results  than  the  traditional  approach. 
Compared  to  the  traditional  approach,  the  new  approach  performs  well,  especially  at  the  surface. 
The  consistency  of  the  new  approach  can  be  seen  in  Figures  42  and  43,  where  the  results  predicted 
by  simulations  follow  the  trends  with  the  experimental  results.  However  for  the  peak  pressure  of 
6.1  GPa,  the  new  approach  underestimates  at  the  surface.  This  could  be  attributed  to  the  number 
of  shots.  The  number  of  shots  for  peak  pressures  of  5.5  and  8.3  GPa  is  only  one,  while  the  number 
of  shots  for  6.1  GPa  is  three.  The  constitutive  material  model  is  not  optimized  for  multiple  shots. 
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Residual  Stress  [MPa] 
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Figure  41:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  JC  Model  with  5.5  Peak  Pressures 
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Figure  42:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  JC  Model  with  6.1  Peak  Pressures 
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Figure  43:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  JC  Model  with  8.3  Peak  Pressures 


The  common  feature  of  all  three  peak  pressures  in  the  new  approach  is  underestimation  of  the 
residual  stresses  predictions  compared  to  the  experimental  results.  This  approach  could  be 
considered  as  a  conservative  design.  The  results  of  the  residual  stress  using  the  ZA  model  for  the 
three  peak  pressures  of  5.5,  6.1,  and  8.3  GPa  are  shown  in  Figures  44,  45,  and  46  respectively. 
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Figure  44:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  ZA  Model  with  5.5  Peak  Pressures 
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Figure  45:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  ZA  Model  with  6.1  Peak  Pressures 
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Figure  46:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  ZA  Model  with  8.3  Peak  Pressures 


The  consistency  of  the  approach  is  evident  in  these  figures  compared  to  experimental  results.  The 
new  approach  performs  well  at  the  surface  for  peak  pressure  of  5.5  GPa  and  6.1  GPa,  while  it 
overestimates  at  the  surface  for  the  peak  pressure  of  8.3  GPa.  On  the  other  hand,  the  traditional 
approach  overestimates  for  all  the  peak  pressures  compared  to  the  experimental  results.  Figures 
47,  48,  and  49  show  the  residual  stress  comparison  between  the  experiment  and  simulations  for 
the  KHL  model  at  peak  pressures  of  5.5,  6.1,  and  8.3  GPa  respectively.  As  seen  in  the  previous 
two  models,  the  KHL  model  performs  well  compared  to  the  experimental  results  for  three 
different  peak  pressures  including  the  surface.  The  same  model  using  the  traditional  approach 
over  predicts  at  the  surface  for  peak  pressures  of  5.5  and  6.1  GPa.  This  shows  that  the 
optimization-based  approach  is  consistent  in  predicting  the  residual  stresses  using  three  different 
models.  In  most  cases  in  this  work,  it  performs  better  than  the  traditional  approach.  However,  we 
cannot  include  this  for  every  case  because  the  residual  stress  prediction  based  on  the  traditional 
approach  is  dependent  on  the  amount  and  quality  of  the  experimental  data  available. 


65 

Approved  for  public  release;  distribution  unlimited. 


Residual  Stress  [MPa] 


100 

0 

-100 

-200 

-300 

-400 

-500 

0  0.5  1  1.5  2 

Depth  [mm]  at  a  distance  of  1  mm  from  center 

Figure  47:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  KHL  Model  with  5.5  Peak  Pressures 


66 

Approved  for  public  release;  distribution  unlimited. 


Residual  Stress  [MPa] 


100 

0 

-100 

-200 

-300 

-400 

-500 

-600 

-700 

0  0.5  1  1.5  2 

Depth  [mm]  at  a  distance  of  1  mm  from  center 

Figure  48:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  KHL  Model  with  6.1  Peak  Pressures 
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Figure  49:  Residual  Stress  Comparison  for  Ti-6A1-4V  for  the  KHL  Model  with  8.3  Peak  Pressures 


Table  18  shows  the  least  square  error  using  Equation  12  and  the  normalized  error  between  the 
experimental  data  and  the  simulation  results  for  all  models.  The  simulation  predictions  include 
the  traditional  approach  (non-optimized)  and  the  optimization-based  approach. 
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This  table  provides  a  quantitative  comparison  between  the  two  approaches.  In  Table  18,  the  least 
squared  error  is  normalized  to  the  least  value  that  occurs  for  the  JC  model  using  the  optimized 
approach  at  the  peak  pressure  of  5.5  GPa.  The  JC  and  the  KHL  models  perform  very  well 
compared  to  the  ZA  model  for  all  the  peak  pressures  for  both  approaches.  It  can  be  seen  that  the 
maximum  error  for  the  optimization-based  approach  is  8.72  times  the  least  error  value  while  it  is 
12.04  times  the  least  error  value  for  the  traditional  approach. 

5.6  Section  Summary  and  Conclusions 

In  this  chapter,  an  inverse  optimization-based  approach  is  used  to  obtain  model  constants  when 
very  little  or  no  experimental  data  of  stress-strain  curves  is  available.  The  optimization-based 
approach  is  shown  to  predict  residual  stresses  that  are  consistent  with  experimental  results.  The 
consistency  of  the  approach  is  shown  by  validating  for  two  materials  of  Inconel®7 1 8  and 
Ti-6A1-4V.  LP  experiments  were  performed  with  a  Nd-glass  laser  for  Inconel®718  at  four 
different  energy  densities  and  the  residual  stress  measurements  were  made  using  an  x-ray 
diffraction  method.  For  the  Inconel®718,  the  JC  and  the  KHL  models  predicted  the  trends  and 
the  simulation  results  are  in  agreement  with  the  experimental  results  for  the  lower  two  peak 
pressures.  The  JC  and  the  KHL  models  are  shown  to  perform  better  than  the  ZA  model  in 
prediction  of  residual  stresses  compared  to  experimental  data  for  Ti-6A1-4V. 
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6.0  PROBABILISTIC  FRAMEWORK  FOR  LASER  PEENED  RESIDUAL  STRESS 
FIELD  BY  PROPAGATION  OF  REGRESSION  UNCERTAINTY 


This  chapter  develops  a  methodology  to  quantify  the  uncertainty  in  residual  stress  field,  i.e., 
residual  stress  value  along  the  depth  induced  by  the  LP  process  simulation.  The  LP  process 
simulation  is  described  in  the  context  of  uncertainty  quantification.  The  cause  of  uncertainty  in 
simulation  of  LP  is  discussed.  The  input  uncertainty  is  the  regression  estimates  obtained  from  the 
non-linear  regression  analysis.  The  effect  of  these  uncertainties  on  the  residual  stress  field  is 
presented. 

6.1  Deterministic  LP  Simulation  Procedure 

The  simulation  procedure  to  predict  residual  stress  field  induced  by  the  LP  process  has  been 
discussed  in  the  previous  chapters.  The  simulation  procedure  is  described  in  this  section  with  the 
purpose  of  making  distinctions  between  deterministic  simulation  procedure  with  stochastic 
simulation  procedure.  A  schematic  of  the  simulation  procedure  is  shown  in  Figure  50. 


Figure  50:  Deterministic  Simulation  Procedure  of  the  LP  Process 


The  deterministic  simulation  procedure  can  be  divided  into  three  components: 

1 .  Obtain  the  estimates  of  material  model  constants  from  a  non-linear  regression  analysis  based 
on  limited  experimental  data 

2.  Input  the  model  constant  estimators  and  other  process  parameters  such  as  peak  pressure,  spatial 
and  temporal  pressure  profile,  and  spot  size  into  FEA 

3.  Compare  the  predicted  residual  stress  field  at  different  depths  with  the  experimental  data 
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The  material  model  constant  estimates  act  as  one  set  of  inputs  to  the  FEA  and  residual  stress  field 
is  obtained  as  the  response  from  the  FEA. 

6.2  Cause  and  Effect  of  Uncertainty  in  LP  Simulation 

During  the  simulation  process,  a  material  model  is  required  which  acts  as  one  of  the  inputs  to 
FEA.  The  residual  stress  field  obtained  is  considered  as  the  response.  A  material  model  is  an 
equation  relating  stresses  to  strains  at  different  strain  rates  and  predicts  the  behavior  of  the 
material  under  different  loading  conditions.  Each  material  model  has  certain  parameters  or  model 
constants  that  must  be  estimated.  These  model  constants  are  determined  using  non-linear 
regression  and  least  squares  approach  is  typically  used  to  estimate  the  model  constants.  An  error 
term  in  the  non-linear  model  is  assumed  to  be  random  and  therefore  estimates  of  the  model 
constants  obtained  are  also  random.  This  randomness  in  the  model  constant  estimates  causes 
uncertainty  in  the  input  and  is  termed  as  regression  uncertainty.  This  uncertainty  represents  the 
variability  of  the  material  model  to  predict  the  material  behavior.  The  uncertainty  in  the  input  is 
propagated  through  the  system  (FEA)  or  black  box  to  obtain  the  uncertainty  in  the  response 
(residual  stress  field).  The  uncertainty  in  the  response  is  expressed  as  confidence  band  on  the 
residual  stress  field. 

6.3  Probabilistic  Framework  for  LP  Simulation 

In  this  research,  a  framework  is  developed  to  address  the  uncertainty  in  the  model  constant 
estimates  that  leads  to  variability  in  the  residual  stress  field.  The  schematic  representation  of  the 
framework  is  shown  in  Figure  5 1 . 
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Figure  51:  Probabilistic  Framework  for  LP  Simulation 


The  framework  is  divided  into  three  components: 

1.  Regression  (Input)  Uncertainty  Quantification 

2.  Uncertainty  Propagation 

3.  Residual  Stress  Field  Confidence  Band 

6.3.1  Regression  (Input)  Uncertainty 

Regression  analysis  is  a  technique  from  statistics  that  can  be  used  to  investigate  the  relationship 
between  a  dependent  variable  and  one  or  more  independent  variables.  Applications  of  the 
regression  technique  are  found  in  many  fields  including  engineering,  physical  sciences, 
economics,  management,  and  biological  sciences.  According  to  Montgomery  and  Peck  [84], 
regression  analysis  is  the  most  widely  used  statistical  technique.  Regression  analysis  can  be  used 
for  purposes  such  as  data  description,  parameter  estimation,  prediction  and  estimation,  and 
control.  Regression  analysis  can  be  categorized  into  two  divisions: 

1 .  Linear  Regression  Analysis 

2.  Non-linear  Regression  Analysis 


73 

Approved  for  public  release;  distribution  unlimited. 


The  simplest  linear  regression  model  is  shown  in  Equation  14. 


y  —  A)  +  Pix + s 


(14) 


where  x  is  the  independent  variable,  /3o  and  /3i  are  the  regression  coefficients  that  need  to  be 
estimated,  y  is  the  dependent  variable,  and  e  is  a  random  error.  Equation  14  is  considered  a  linear 
regression  model  because  the  equation  is  linear  in  coefficients.  The  least  squares  method  is 
typically  used  to  estimate  the  regression  coefficients.  The  error  term  is  assumed  to  be  Nll)(0.  a2), 
normal  and  independently  distributed  with  zero  mean  and  constant  variance  of  o2.  The 
assumption  of  zero  mean  in  the  error  is  primarily  used  for  simplifying  the  number  of  unknown 
regression  coefficients.  A  non-zero  mean  is  absorbed  in  the  intercept  Q3o)  of  the  model.  The 
assumption  of  constant  variance  of  value  a2  simply  states  that  the  observed  experimental  data 
values  obtained  are  equally  unreliable.  Further  analysis  can  be  performed  to  test  this  assumption 
and  variance-stabilization  transformation  can  also  be  implemented  [85].  The  least  squares 
estimator  of  vector  /3  is  shown  in  Equation  15: 

&  =  (XTX)-lXTy  (15) 


where  X  is  a  matrix  of  levels  of  the  independent  variables  and  a  column  of  ones  for  the  intercept 
and  y  is  the  vector  of  observations.  The  size  of  y  is  n  x  1,  where  n  is  the  number  of  observations. 
The  size  of  X  is  n  x  k,  where  k  represents  the  number  of  regression  coefficients  which  includes  the 
intercept  term  /3q.  The  size  of  the  vector  j8  is  k  x  1.  The  fitted  regression  model  is  shown  in 
Equation  16: 

y  =  Xp  (16) 

where  y  is  known  as  fitted  values  or  predicted  values.  The  difference  between  actual  observation 
yi  and  fitted  value  y,  is  the  residual,  et  —  y,  —  y,.  The  vector  of  residuals  n  x  1  is  e  =  y  —  y.  The 
formula  for  estimating  the  error  variance  (c2)  is  shown  in  Equation  17: 


xTi 

°  n-k 


(17) 


The  variance  property  of  j8  is  expressed  as  a  covariance  matrix  obtained  by 

Cov(j 8 )  =  o2  (XTX )  1  [86].  Error  variance  estimated  in  Equation  17  is  used  for  the  error  variance 
term  to  obtain  the  estimate  of  the  covariance  matrix.  The  size  of  the  covariance  matrix  is  k  x  k. 
The  covariance  matrix  is  a  symmetric  matrix  whose  ith  diagonal  element  is  the  variance  of  the 
individual  regression  coefficient  estimate  j8,  and  i jth  element  is  the  covariance  between  estimates 
j8,  and  j3j.  Assuming  a  normal  error,  the  regression  coefficient  estimators  /3  will  have  multivariate 
normality  distribution.  Confidence  intervals  for  individual  regression  coefficients  can  be  found  in 
Montgomery’s  work  [86]. 


Engineers  experience  non-linear  regression  models  in  day  to  day  life.  The  common  notation  for  a 
non-linear  regression  model  is  shown  in  Equation  18: 

y  =  g(x;po,Pi,---,fik)  +  e  (!8) 
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where  x  is  the  vector  of  independent  variables  and  /3,  i  =  1.2...../:  are  the  regression  coefficients 
that  need  to  be  estimated.  The  term  ’non-linear’  is  used  because  the  function  g  is  non-linear  in  the 
coefficients  (j3;).  Unlike  linear  regression  models,  there  is  no  closed  form  solution  available  to 
estimate  the  regression  coefficients  in  a  non-linear  model.  Instead,  iterative  optimization 
techniques  such  as  Gauss-Newton  algorithm  is  used.  The  basic  idea  of  the  Gauss-Newton 
algorithm  is  approximating  the  function  g  by  a  first  order  Taylor  series  expansion.  This  linear 
approximation  is  used  in  least  squares  approach.  A  Jacobian  matrix,  which  is  composed  of  partial 
differential  of  function  g  with  respect  to  each  parameter  /3;,  is  obtained.  An  initial  value  is 
proposed  for  the  regression  coefficients  and  the  coefficients  are  updated  until  convergence  is 
reached.  The  details  of  the  algorithm  can  be  found  in  the  work  of  Ratkowsky  [87]. 

The  covariance  matrix  [iff]  between  the  regression  coefficient  estimates  can  be  estimated  from  the 
linearized  approximation  of  the  function  g.  The  partial  differential  of  function  g  with  respect  to 
each  regression  estimate  /3 j  for  an  independent  variable  is  denoted  by  Equation  19: 


and  [C]  =  {cij}-  Compared  to  the  linear  regression  analysis,  matrix  [C]  is  analogous  to  matrix 
[A].  The  size  of  the  matrix  [C]  is  n  x  k,  where  n  is  number  of  observations  and  k  represents  the 
number  of  regression  coefficients.  The  equation  for  the  estimated  covariance  between  the 
regression  estimates  is  shown  in  Equation  20: 

[y}  =  62[CTC\~l  (20) 

where  <7  2  is  the  estimated  error  variance  as  shown  in  Equation  17. 

The  estimates  of  the  regression  coefficients  are  approximately  normal.  The  approximate  normal 
distribution  gets  closer  to  the  asymptotic  limit  as  the  sample  size  gets  larger  based  on  law  of  large 
numbers  and  central  limit  theorem.  For  a  limited  sample  size,  which  is  mostly  the  case  in 
engineering,  it  is  desirable  to  validate  the  normality  assumption.  In  this  research  an  approach 
from  statistics  known  as  “bootstrap  for  regression”  is  used  to  verify  the  normality  assumptions  of 
the  regression  coefficient  estimates  [88].  The  following  section  details  bootstrapping. 

6.3.2  What  is  Bootstrapping? 

Bootstrapping  is  a  modern,  computer-intensive  statistical  technique.  It  has  been  initiated  by  Efron 
and  the  first  monograph  was  published  in  1993  [89].  Bootstrapping  can  be  used  to  estimate  the 
standard  errors  of  means,  variances  as  well  as  more  complicated  statistics.  It  can  also  be  used  to 
construct  hypothesis  tests.  Bootstrapping  falls  into  the  category  of  sampling  with  replacement 
(except  few  exceptions).  Bootstrapping  provides  additional  details  than  a  simple  statistic  from  a 
sample.  Figure  52  shows  a  generic  flow  chart  of  the  bootstrapping  method.  The  figure  can  be 
easily  understood  with  an  example.  Suppose  the  initial  sample  data  of  an  observed  quantity  is  10, 
13,  15,  18,  and  20.  ’Mean’  of  the  data  can  be  obtained  and  represents  the  sample  statistic.  The 
mean  of  the  above  data  is  15.2.  Bootstrapping  can  provide  us  additional  information  such  as  the 
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Figure  52:  Schematic  Representation  of  Bootstrap  Method 

distribution  of  the  mean.  It  can  be  assumed  that  there  are  five  balls  in  a  basket  numbered  10,  13, 
15,  18,  and  20.  A  ball  is  randomly  sampled  and  then  placed  back  in  the  basket.  This  is  known  as 
sampling  with  replacement.  There  are  again  all  five  balls  in  the  basket  when  the  next  ball  is 
randomly  sampled.  The  procedure  is  repeated  five  times  to  obtain  the  first  bootstrap  sample.  The 
first  bootstrap  sample  obtained  using  MATLAB  [90]  is  18,  13,  20,  13,  and  18  and  the  mean  of  the 
bootstrap  sample  is  16.4.  As  it  can  be  seen  from  the  first  bootstrap  sample,  the  numbers  13  and  18 
appears  twice,  the  numbers  10  and  15  does  not  appear  at  all.  Another  possible  bootstrap  sample 
can  be  that  all  five  samples  can  be  the  same  number.  In  this  manner  p  bootstrap  samples  are 
generated  from  the  initial  sample.  The  typical  values  of  p  range  from  500  to  5000.  Figure  53 
shows  the  histogram  of  the  3000  bootstrap  samples.  The  number  of  bootstrap  samples  (3000)  is 
chosen  large  enough  so  that  the  results  are  invariant  of  the  number  of  bootstrap  samples.  The 
variation  of  the  mean  that  is  centered  around  the  sample  mean  can  be  seen  from  the  histogram 
data.  The  power  of  the  method  lies  in  the  fact  that  additional  samples  are  being  generated  from 
the  initial  sample.  This  indicates  that  no  new  experiments  are  conducted  to  obtain  the  additional 
information.  In  the  above  example,  the  sample  statistic  was  the  mean  of  the  data,  but  the 
technique  can  also  be  used  for  other  sample  statistics.  Bootstrap  method  has  been  used  by 
Picheny  et.  al.  to  conservatively  estimate  the  reliability  of  engineering  structures.  Probability  of 
failure  of  a  composite  panel  is  estimated  using  different  methods  and  the  results  are  compared 
with  the  bootstrap  method  [91]. 

6.3.2. 1  Bootstrapping  for  Regression 

Bootstrapping  can  also  be  used  for  regression.  There  are  three  ways  to  bootstrap  a  regression 
analysis  [92],  which  are  divided  into  two  categories:  (i)  parametric  bootstrap  and  (ii) 
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Figure  53:  Histogram  of  Mean  of  the  Sample  Data 


non-parametric  bootstrap.  Example  10-1  from  Montgomery’s  book  [86]  in  Chapter  10  is  used  to 
demonstrate  the  bootstrap  methods.  Sixteen  observations  of  an  experiment  on  the  viscosity  of  a 
polymer  (y)  and  two  process  variables,  reaction  temperature  (x\ )  and  catalyst  feed  rate  (*2),  are 
shown  in  Table  19.  A  linear  regression  model  as  shown  in  Equation  21  is  used  to  fit  the 
experimental  data. 


y  —  A)  +  +  P2X2  +  £ 


(21) 
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Table  19:  Viscosity  Data  of  Demonstration  Example  [86] 


Observation 

Temperature  (jci) 

Catalyst  Feed  Rate  (*2) 

Viscosity  (y) 

1 

80 

8 

2256 

2 

93 

9 

2340 

3 

100 

10 

2426 

4 

82 

12 

2293 

5 

90 

11 

2330 

6 

99 

8 

2368 

7 

81 

8 

2250 

8 

96 

10 

2409 

9 

94 

12 

2364 

10 

93 

11 

2379 

11 

97 

13 

2440 

12 

95 

11 

2364 

13 

100 

8 

2404 

14 

85 

12 

2317 

15 

86 

9 

2309 

16 

87 

12 

2328 

The  least  squares  estimates  of  /3  obtained  using  Equation  15  are  shown  in  Table  20: 


Table  20:  Regression  Coefficients  for  Linear  Regression  Model 


Parameter 

Value 

A) 

1566.08 

A 

7.62 

A 

8.58 

The  variance  estimated  for  the  error  (cr2)  using  the  Equation  17  is  16.36.  Table  21  shows  the 
observations,  the  predicted  values,  and  the  corresponding  residuals  or  errors.  The  predicted  values 
are  calculated  using  Equation  16. 
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Table  21:  Observations,  Predicted  Values,  and  Residuals  for  Demonstration  Example  [86] 


# 

Observation  (>’,•) 

Predicted  value  (y,) 

Residual  (<?,•  =  yi  -  yi) 

1 

2256 

2244.5 

11.5 

2 

2340 

2352.1 

-12.1 

3 

2426 

2414.1 

11.9 

4 

2293 

2294.0 

-1.0 

5 

2330 

2346.4 

-16.4 

6 

2368 

2389.3 

-21.3 

7 

2250 

2252.1 

-2.1 

8 

2409 

2383.6 

25.4 

9 

2364 

2385.5 

-21.5 

10 

2379 

2369.3 

9.7 

11 

2440 

2416.9 

23.1 

12 

2364 

2384.5 

-20.5 

13 

2404 

2396.9 

7.1 

14 

2317 

2316.9 

0.1 

15 

2309 

2298.8 

10.2 

16 

2328 

2332.1 

-4.1 

Parametric  Bootstrap 

The  parametric  bootstrap  is  based  on  the  belief  that  the  error  term  in  the  regression  analysis 
shown  in  Equation  14  or  18  is  Nll)(0.  a2),  normal  and  independently  distributed  with  zero  mean 
and  constant  variance.  The  estimates  of  the  regression  coefficients  j  are  obtained  from  the 
regression  analysis  as  shown  in  Table  20.  The  variance  of  the  error  is  estimated  using  Equation  17 
and  the  value  obtained  is  16.36  for  the  above  demonstration  example.  It  is  termed  as  a2.  Samples 
of  the  error  are  generated  based  on  normal  distribution  £*  ~  A  (0,  cf2) .  The  method  is  termed  as 
’parametric’  because  the  bootstrap  samples  are  generated  from  a  distribution  (e.g.  normal)  instead 
of  resampling  from  the  original  data.  The  number  of  samples  generated  is  equal  to  the  number  of 
observations.  In  the  above  example,  since  16  observations  are  available,  the  number  of  normal 
random  samples  of  error  is  also  16.  New  observation  is  generated  by  adding  the  sampled  error  to 
the  predicted  values,  y*  —  y,-  +  £*,  i  =  1,2, ...  ,n,  where  n  is  the  number  of  observations.  In  this 
way  a  new  set  of  data  (x,-,y*)  is  generated.  Using  the  newly  generated  data,  regression  analysis  is 
performed  to  obtain  the  estimates  of  regression  coefficients.  The  procedure  is  repeated  p  times, 
where  p  is  the  number  of  bootstraps.  The  method  is  not  valid  if  the  error  is  not  normal  and  this 
can  be  a  disadvantage  of  the  parametric  bootstrap  approach.  Step  by  step  procedure  of  the  method 
for  the  viscosity  of  polymer  example  discussed  in  the  previous  section  is  shown  in  the  following. 

Random  samples  of  a  normal  distribution  for  the  error  with  a  mean  of  zero  and  variance  of  16.36 
(estimated  variance)  is  generated.  The  size  of  the  sample  is  equal  to  number  of  observations, 
which  is  16.  The  errors  are  added  to  the  fit  or  predicted  values  to  generate  a  new  observation. 
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Table  22  shows  one  example  of  the  first  randomly  sampled  error,  the  fit,  and  the  generated 
observation. 


Table  22:  Randomly  Sampled  Error,  Predicted  Values,  and  Generated  Observations 


# 

Random  sample  of  error  (£■  ) 

Predicted  value  (y,-) 

Generated  observation  (y*  =  y,  4-  e* ) 

1 

2.0667 

2244.5 

2246.57 

2 

-43.5123 

2352.1 

2308.59 

3 

5.3382 

2414.1 

2419.44 

4 

-8.8267 

2294.0 

2285.17 

5 

20.0801 

2346.4 

2366.48 

6 

14.4624 

2389.3 

2403.93 

7 

8.7680 

2252.1 

2260.87 

8 

12.9395 

2383.6 

2396.54 

9 

24.0032 

2385.5 

2409.50 

10 

-21.8989 

2369.3 

2347.40 

11 

1.5361 

2416.9 

2418.44 

12 

-7.9748 

2384.5 

2376.53 

13 

-23.4989 

2396.9 

2373.40 

14 

-15.8406 

2316.9 

2301.06 

15 

19.0240 

2298.8 

2317.82 

16 

14.7254 

2332.1 

2346.83 

The  first  bootstrap  estimates  of  /3  obtained  for  the  newly  generated  observation  using  Equation  15 

are  /3q  =  1579.12,  jSj*  =  7.37,  and  /3|  =  9.56.  The  procedure  of  generating  normal  random 
samples  of  error  is  repeated,  and  for  every  sample  a  bootstrap  estimate  of  the  regression 
coefficients  is  obtained.  Figures  54,  55,  and  56  show  the  histogram  of  3000  bootstrap  estimates  of 
j8q  ,  /3|\  and  /3|  respectively.  The  number  of  bootstrap  samples  (3000)  is  chosen  large  enough  so 
that  the  results  are  invariant  of  the  number  of  bootstrap  samples.  Additional  information  about  the 
estimates  is  obtained  at  no  further  cost  of  performing  an  actual  experiment. 

Non-parametric  Bootstrap 

There  are  two  methods  in  non-parametric  bootstrap.  They  are  termed  as  (i)  observation 
resampling  or  random-x  resampling  and  (ii)  residual  resampling  or  fixed-x  resampling.  Each 
method  is  described  and  comparison  between  the  two  are  made  later.  The  method  is  called  as 
’non-parametric’  because  the  bootstrap  samples  are  not  generated  from  a  known  distribution, 
which  is  the  case  for  a  parametric  bootstrap.  Instead,  the  samples  are  generated  by  resampling  the 
original  data. 

Random- x  Bootstrap:  Demonstration  example  discussed  in  the  previous  section  is  used  to 
explain  the  method.  The  viscosity  data  of  a  polymer  is  a  vector  of  observations  (y)  and  the 
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Figure  54:  Parametric  bootstrap  histogram  of  j8q  for  the  demonstration  example 
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Figure  55:  Parametric  bootstrap  histogram  of  jSffor  the  Demonstration  Example 
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Figure  56:  Parametric  Bootstrap  Histogram  of  j8|  for  the  Demonstration  Example 


model-matrix  (X)  can  be  written  as  shown  in  Equation  22: 


"  1 

80 

8 

'  2256  ' 

1 

93 

9 

2340 

1 

100 

10 

2426 

1 

82 

12 

2293 

1 

90 

11 

2330 

1 

99 

8 

2368 

1 

81 

8 

2250 

1 

96 

10 

2409 

1 

94 

12 

y  = 

2364 

1 

93 

11 

2379 

1 

97 

13 

2440 

1 

95 

11 

2364 

1 

100 

8 

2404 

1 

85 

12 

2317 

1 

86 

9 

2309 

1 

87 

12 

16x3 

2328 

A  matrix  Z  is  created  by  combining  the  model-matrix  X  and  vector  of  observations  y.  The  first 
column  of  the  Z  matrix  is  the  vector  of  observations  and  the  remaining  columns  are  from  the 
model-matrix  X.  The  Z  matrix  obtained  for  the  demonstration  example  from  the  Equation  22  is 
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shown  in  Equation  23: 


Z  = 


2256 

1 

80 

8 

2340 

1 

93 

9 

2426 

1 

100 

10 

2293 

1 

82 

12 

2330 

1 

90 

11 

2368 

1 

99 

8 

2250 

1 

81 

8 

2409 

1 

96 

10 

2364 

1 

94 

12 

2379 

1 

93 

11 

2440 

1 

97 

13 

2364 

1 

95 

11 

2404 

1 

100 

8 

2317 

1 

85 

12 

2309 

1 

86 

9 

2328 

1 

87 

12 

(23) 


16x4 


Each  row  of  the  Z  matrix  is  considered  as  a  sample  data.  Bootstrap  samples  from  the  rows  of  the 
Z  matrix  are  generated  and  is  termed  as  Z* .  Equation  24  shows  an  example  of  a  Z*  matrix. 


Z*  = 


2340 

1 

93 

9 

2250 

1 

81 

8 

2250 

1 

81 

8 

2309 

1 

86 

9 

2440 

1 

97 

13 

2317 

1 

85 

12 
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1 

90 

11 

2404 

1 

100 

8 

2426 

1 

100 

10 

2309 

1 

86 

9 

2379 

1 

93 

11 

2328 

1 

87 

12 

2364 

1 

95 

11 

2340 

1 

93 

9 

2250 

1 

81 

8 

2330 

1 

90 

11 

(24) 


J  16x4 

The  first  row  the  Z*  matrix  is  obtained  from  randomly  sampling  from  one  of  the  rows  of  the  Z 
matrix.  It  can  be  seen  from  Equation  24,  the  first  row  of  the  Z*  matrix  is  the  second  row  of  the  Z 
matrix.  The  sampled  row  from  the  Z  matrix  is  replaced  back  so  that  all  the  rows  of  the  Z  matrix 
are  available  to  randomly  sample  for  the  second  row  for  the  Z*  matrix.  In  this  manner  all  the 
sixteen  rows  of  the  Z*  matrix  are  generated.  It  can  be  seen  from  the  comparison  between 
Equations  23  and  24  that  rows  1,  4,  6,  8,  and  9  of  Z  matrix  don’t  appear  in  the  Z*  matrix  at  all, 
where  as  row  7  appears  thrice,  rows  5  and  15  appears  twice,  and  the  remaining  rows  appear  once. 
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The  matrix  Z*  can  be  decomposed  into  bootstrap  model  matrix  X*  and  bootstrap  observation 
vector  y*  as  shown  in  Equation  25: 
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9 
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81 
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1 

81 
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86 

9 
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85 

12 

1 

90 

11 
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100 
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1 

86 
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1 

93 

11 

1 

87 
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1 

95 

11 

1 

93 

9 

1 

81 

8 

1 

90 

11 

'  2340  " 
2250 
2250 
2309 
2440 
2317 
2330 
2404 
2426 
2309 
2379 
2328 
2364 
2340 
2250 


(25) 


Bootstrap  estimates  of  the  coefficients  obtained  using  Equation  15  are  /3q  =  1536.85,  /3f  =  7.926, 
and  j8|  =  8.671.  The  procedure  is  repeated  3000  times  and  each  time  the  bootstrap  estimates  are 
saved.  The  number  of  bootstrap  samples  (3000)  is  chosen  large  enough  so  that  the  results  are 
invariant  of  the  number  of  bootstrap  samples.  Figures  57,  58,  and  59  show  the  histograms  of  the 
bootstrap  estimates  of  /3()  .  /3j\  and  /3|  respectively. 

Fixed-x  Bootstrap:  Regression  analysis  is  performed  using  Equation  15  for  a  linear  model  to 
obtain  the  estimates.  Table  20  shows  the  obtained  estimates  for  the  demonstration  example.  Table 
21  shows  the  observations  (y,-)  ,i  —  1, 2, ...  n,  the  predicted  values  or  the  fit  (y,-)  ,i  —  1, 2, ...  n,  and 
the  errors  (c,-) ,  i  —  1,2, ...  n.  n  is  the  number  of  observations.  A  bootstrap  sample  from  the 
obtained  errors  (e,-)  is  generated  and  is  termed  as  e*,i=  1,2 The  number  of  bootstrap 
errors  is  equal  to  the  number  of  observations.  New  observations  are  generated  by  adding  the  fit  to 
the  bootstrap  errors,  y*  =  Si  +  e*,  i  —  1,2,...  n.  In  this  way  a  new  data  set  is  generated.  Using  the 
newly  generated  data  set,  regression  analysis  is  performed  to  obtain  bootstrap  estimates  of  the 
coefficients.  The  procedure  of  bootstrapping  the  errors  is  repeated  and  regression  analysis  is 
performed  to  obtain  the  bootstrap  estimates.  The  bootstrap  estimates  are  stored  for  each 
regression  analysis  to  obtain  a  distribution  of  the  regression  estimates.  The  difference  between  the 
parametric  bootstrap  and  the  fixed-x  bootstrap  is  that,  for  a  parametric  bootstrap,  the  errors  are 
sampled  from  a  normal  distribution  of  zero  mean  and  62  variance,  while  the  errors  are  sampled 
with  replacement  from  the  original  regression  for  the  fixed-x  bootstrap.  Figure  60  shows  the  flow 
chart  of  the  fixed-x  resampling  procedure.  Step  by  step  procedure  of  the  method  for  the 
demonstration  example  is  shown  in  the  following. 

Table  21  shows  the  observations,  predicted  values,  and  the  errors  of  the  16  data  points  from  the 
original  regression  analysis.  The  errors  from  the  original  regression  analysis  are  sampled  with 
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Figure  57:  Random-x  Bootstrap  Histogram  of  j8q  for  the  Demonstration  Example 
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Figure  58:  Random-x  Bootstrap  Histogram  of  /3*  for  the  Demonstration  Example 
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Figure  59:  Random-x  Bootstrap  Histogram  of  /3|  for  the  Demonstration  Example 


Figure  60:  Fixed-x  Bootstrap  Procedure  Flowchart 
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replacement.  The  size  of  the  error  vector  is  same  as  the  number  of  observation,  which  is  16  in  the 
demonstration  example  Table  23  shows  an  example  of  the  bootstrap  errors,  fit  from  the  original 
regression  analysis,  and  the  generated  observation.  The  generated  observation  is  obtained  by 
adding  the  fit  to  the  bootstrap  errors.  Comparing  the  bootstrap  error  column  in  Table  23  to  the 
error  column  in  Table  21,  it  can  be  seen  that  the  errors  11.9,  -2.1,  7.1,  0.1,  and  10.2  from  Table  21 
did  not  appear  in  the  bootstrap  errors,  while  the  errors,  -1.0  and  -21.5  appeared  three  times  and 
the  error  23.1  appeared  twice. 


Table  23:  Bootstrap  Error,  Predicted  Values,  and  Generated  Observations  for  Fixed-x  Bootstrap 
Method 


# 

Bootstrap  error  (e*) 

Predicted  value  (y;) 

Generated  observation  ( y *  —  y,-  +  £;* ) 

1 

-12.1 

2244.5 

2232.4 

2 

-1.0 

2352.1 

2351.1 

3 

11.5 

2414.1 

2425.6 

4 

-20.5 

2294.0 

2273.5 

5 

25.4 

2346.4 

2371.8 

6 

-4.1 

2389.3 

2385.2 

7 

23.1 

2252.1 

2275.2 

8 

-21.3 

2383.6 

2362.3 

9 

-1.0 

2385.5 

2384.5 

10 

-16.4 

2369.3 

2352.9 

11 

-21.5 

2416.9 

2395.4 

12 

-1.0 

2384.5 

2383.5 

13 

-21.5 

2396.9 

2375.4 

14 

23.1 

2316.9 

2340 

15 

9.7 

2298.8 

2308.5 

16 

-21.5 

2332.1 

2310.6 

The  bootstrap  estimates  obtained  from  the  newly  generated  data  set  using  Equation  15  are 
j8q  =  1633.34,  j8j*  =  7.0056,  and  /3|  =  7.19.  The  procedure  of  generating  bootstrap  errors  is 
repeated,  and  for  every  sample  of  the  bootstrap  estimate  of  the  regression  coefficients  is  obtained 
and  saved.  The  procedure  is  repeated  3000  times  to  obtain  the  histograms  of  the  bootstrap 
estimates.  Figures  61,  62,  and  63  show  the  bootstrap  histograms  of  j8q  ,  /3j\  and  j89*  respectively. 

The  fixed-x  resampling  technique  is  used  when  the  model-matrix  [X]  is  fixed.  Efron  [89]  claims 
that  the  two  approaches  are  asymptotically  equivalent.  The  random-x  method  is  less  sensitive  to 
assumptions  concerning  independence  or  error  term.  In  this  work,  the  fixed-x  method  is  suitable 
because  the  model-matrix  [X]  is  fixed. 

At  the  end  of  the  input  uncertainty  quantification  step,  the  covariance  matrix  [y/]  between  the 
regression  estimates  which  represents  the  input  uncertainty  is  obtained  along  with  the  estimates 
itself.  The  size  of  the  covariance  matrix  is  k  x  k,  where  k  is  the  number  of  regression  estimates. 
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Figure  61:  Fixed-x  Bootstrap  Histogram  of  j8q  for  the  Demonstration  Example 
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Figure  62:  Fixed-x  Bootstrap  Histogram  of  j8j*  for  the  Demonstration  Example 
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Figure  63:  Fixed-x  Bootstrap  Histogram  of  /32  for  the  Demonstration  Example 
6.3.3  Uncertainty  Propagation 

Propagation  of  the  uncertainty  is  needed  to  account  its  effect  on  the  response.  In  a 
simulation-based  environment  such  as  FEA,  polynomial  chaos  expansion  or  Taylor  series 
expansion  are  used  to  approximate  the  response.  This  is  performed  to  reduce  the  computational 
cost  because  there  is  no  explicit  relationship  between  the  input  and  the  response.  A  closed  form 
equation  is  used  to  demonstrate  the  derivation  procedure  for  the  input  uncertainty  propagation 
before  the  black  box  (FEA)  propagation  is  presented. 

Closed  Form  Propagation 

Consider  a  simple  function  /,  that  is  a  function  of  two  parameters  j6 o  and  /3i.  Assume 

/  a\  a  a  /V  /V 

y  —  f[t\l 3  )  =  /3()  +  fi\l,  where  l  is  the  independent  variable  and  /3o  and  /3i  are  unbiased  estimates 

of  true  values.  In  analogy  with  the  LP  simulation,  y  can  be  considered  as  residual  stress,  j8o  and 
j8i  are  model  constant  estimates  and  t  is  depth.  The  covariance  function  of  response  at  two 
different  depths  t\  and  t2  is  derived: 


Cov[y(ti),y(i2)]  =  Cov  (po  +  Pih) ,  (j80  +  j8 {t2 


:  E  {  (/3o  +  /3pij  —  (Po  +  Pih)  (jk)  +  Pit2^  —  (Po  +  Pih) 


,  assuming  E  \  /3q  \  —  fio  and 
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E  {  (A)  ^  A))  —  (j$l  —  /3l)  t\  (^Po  -  Po'j  -  (jil  —Pl^t2 


=  E 


/3o-/3o 


+  E 


(/3o —  A>)  (j$i  —  A j  (?i+?2)+£’ 


A~A 


(rite) 


=  600  +  601  (ri  +  te)  +  6n  (rite) 

where  600  is  the  estimate  of  variance  of  /3q,  61 1  is  the  estimate  of  variance  of  /3 1  and  601  is  the 
estimate  of  covariance  between  j8o  and  j}\.  E  (•)  stands  for  expected  value.  The  above  equation 
can  be  written  in  matrix  form: 


y  (*i) 

Var\y(h)\ 

Cov\y(ti),y(t2)] 

.  y(* 2) . 

Cov\y(t2)  ,y  (ti)] 

Var\y(t2 )] 

where  Var  stands  for  variance  of  a  variable  and  Cov  stands  for  covariance  of  a  variable. 
Substituting  for  the  covariance  as  derived  above  and  using  the  property  of  variance  as  self 
covariance: 


y(h) 

obo  +  2/ydoi  +t\&u  600  +  001  (ri  +  te)  +  rin  (*1*2) 

.  yfo) . 

600  +  601  (ri  +  te)  +  6n  (*1*2)  <7oo  +  2riribi  +  t|6n 

(26) 


The  matrix  shown  in  Equation  26  is  decomposed  into  three  matrices  as:  gradient  matrix,  input 
covariance  matrix,  and  transpose  of  gradient  matrix  as  follows: 


Cov 


y(h) 

600  +  ri6oi 

601  +  ri6 11 

'  1  1 

y(* 2) 

dd  +  t2On 

601  +  te6n 

h  te 

\  y(ti) 

1  ri 

- 1 

O 

^01  ’ 

[1  il 

1 - 

to 

1  te 

L  ^01 

<7ll 

1 - 

to* 

1 _ 

(27) 

(28) 


The  matrix  shown  in  Equation  28  can  be  written  as  a  gradient  matrix,  input  covariance  matrix, 
and  transpose  of  the  gradient  matrix  as  shown  in  Equation  29: 


Cov 


y(h) 

y(h) 


df(h)  df{t  1) 

r 

< 

< 

L 

dfih)  df{h) 

dp0  dPi 

Coo  O01 

dp0  dp0 

dfXh)  df{t2) 

obi  on  _ 

dfih)  df{t2) 

.  dp0  rjfi\  _ 

L  dp,  dp,  J 

Cov 


y(ti) 

y{h) 


[F]  M  [ft] 


(29) 


(30) 


where  [F]  is  a  matrix  of  partial  derivatives  of  parameter  estimates  at  two  different  depths  t\  and  h, 
and  [y]  is  the  covariance  matrix  of  the  estimates  of  regression  coefficients  obtained  from  the 


non-linear  regression  analysis.  The  term 


dJM 

L  dpo 


is  the  gradient  of  /  with  respect  to  /3o  at  a  depth 


1\ .  The  covariance  matrix  of  the  response  is  obtained  by  propagating  the  input  uncertainty.  The 
derivation  shown  in  this  section  is  valid  for  a  linear  relationship  between  input  and  output. 
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Black  Box  Propagation 

In  the  above  demonstration  example,  the  relationship  between  random  input  parameters  j8  and  the 
response  y  is  known.  In  a  simulation-based  environment,  such  as  FEA,  the  function  /  is  a  black 
box.  First  order  Taylor  series  expansion  is  used  to  approximate  the  response.  At  any  depth  t,  the 
residual  stress  can  be  approximated  as  shown  in  Equation  31: 


(31) 


where  p  represents  the  estimates  of  material  model  constants  p. 


dm 

dp 


P=  P 


is  a  vector 


representing  the  gradient  of  residual  stress  w.  r.  t.  all  the  regression  coefficients  at  any  depth  t 
evaluated  at  the  regression  estimates  obtained.  The  definition  of  covariance  of  residual  stresses 
between  two  different  depths  t\  and  ti  is  shown  in  Equation  32: 


Cov 


f  (ri  • /^)  -,f  (ft; /^)  =£'|  —  f(tuP^  /  —  /  (?2>/^)  }  (32) 


Substituting  the  Taylor  series  approximation  assumed  in  Equation  31  in  the  Equation  32,  we 
obtain: 


df(t\) 

dP 


P=P 


P~P 


df(ti) 

dp 


P=P 


P-P 


Rearranging  the  terms  in  Equation  33,  we  obtain: 

T 


dm 

dp 


p=p 


p-p  (p-p 


dm 

dp 


p=p 


(33) 


(34) 


The  gradient  matrix  and  its  transpose  shown  in  Equation  34  are  not  random  but  fixed  numbers. 
Therefore  Equation  34  is  simplified  as  shown  in  Equation  35: 


' dm 

T  E{ 

1 

^>j 

1 

Xo>| 

dP 

P=P  1 

LV-  -/  V-  -/  j 

df(t 2) 
dP 


P=P 


(35) 


The  expected  value  of  the  random  term  shown  in  Equation  35  is  by  definition  the  covariance 
matrix  of  the  regression  estimates  that  is  obtained  from  the  regression  analysis.  It  is  defined  as 
shown  in  30  as  [y/j.  Substituting  [y/]  in  Equation  35,  we  obtain: 


df(h ) 

dp 


M 


p=p 


dm 

dp 


p=p 


(36) 
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T 

is  the  partial  derivative  of  residual  stress  at  depth  with  respect  to  estimates  of 

P=P 

model  constants.  For  the  sake  of  simplicity  the  subscript  /3  =  /3  is  not  included  in  rest  of  the 
document.  [\jf]  is  a  symmetric  matrix  with  variance  on  diagonal  terms  and  covariance  on  off 
diagonal  terms.  The  number  of  discrete  depths  that  are  considered  depends  on  the  finite  element 
(FE)  model  mesh.  Let  m  be  the  number  of  nodes  in  the  FE  model.  The  covariance  of  residual 
stress  field  at  each  depth  is  obtained  as  follows: 


dp 


'  f 

(h;P)  ' 

Cov 

f 

= 

C 

./( 

!. 

c 

/  (h ,  /3 ) 

c 

f  ’  (?2’ 

£)] 

•  C 

f{t  i,p) 

1  (  tm  ?  P  ) 

t2'P)  ’ 

■S\ 

V 

f  (t2,Pj 

•  C 

fy*2,i l) 

■>  (  P  ) 

(  -A,  \  /  A  \  /  A,  \  /  -A  \ 

C  J  \  tm,p  \  >  p2>P) 


where  V  [•]  and  C  [•]  represent  the  variance  and  covariance  of  a  variable  respectively.  The  size  of 
the  matrix  is  m  x  m,  with  variance  of  residual  stress  at  each  depth  on  diagonal  terms  and 
covariance  between  depths  on  the  off  diagonal  terms.  Substituting  the  covariance  between 
between  two  depths  that  has  been  derived  before,  we  obtain: 


\df(ti)] 

_  dp  \ 

T 

¥ 

_  dp  \ 

\df{h)  1 
_  dp  \ 

T 

¥ 

\df(t2)] 
[  \ 

\dm] 

.  w  \ 

T 

¥ 

W  . 

r ^ 

p/te) 

[  W  \ 

T 

¥ 

p/(p) 

L  w  \ 

\df(t2 ) 

[  ^  \ 

T 

¥ 

[  \ 

p/te) 

[  w  \ 

T 

¥ 

dfifrn) 

[  dp  \ 

d  f{tm) 

dp  _ 

T 

¥ 

\df(hY 

[  dp  _ 

d  f{tm) 

dP 

T 

¥ 

Seal 

i _ i 

d  f{tm) 

W  . 

T 

¥ 

[  dp  _ 

(37) 


The  above  matrix  can  be  simplified  by  separating  the  common  terms  and  the  resultant  matrices 
obtained  are: 


mi) 

dp 

dJM 

dp 


T 


1  xk 
T 


1  xk 


Mm 


df(t  i)l 

df(tm) 

dp_  \ 

kx  1 

[  dp  \ 

kx  1 

[  dp  \ 

kx  1  . 

[F r]  M  [F] 


df(tm) 

dp 


T 

1  xk 


m  x  k 


where  [F]  is  the  sensitivity  matrix  of  the  residual  stress  with  respect  to  estimates  of  material 
model  constants.  Each  column  of  [F]  is  the  gradient  of  residual  stress  with  respect  to  a  material 
model  constant  estimate  along  the  depth.  The  obtained  covariance  matrix  of  residual  stress 
represents  the  uncertainty  in  the  residual  stress  field  due  to  uncertainty  in  the  model  constant 
estimator.  Central  difference  scheme  is  used  to  determine  the  gradients  of  residual  stress  at  each 
depth  with  respect  to  estimates  of  regression  coefficients  or  model  constants.  The  step  size  for 
each  gradient  is  obtained  after  a  convergence  study.  The  number  of  additional  FE  simulations 
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Figure  64:  Confidence  Band  Validation  Procedure 

needed  using  the  central  difference  scheme  depends  on  the  number  of  regression  coefficient 
estimates.  For  k  estimates,  2k  additional  FE  simulations  are  required. 

6.3.4  Residual  Stress  Field  Confidence  Band 

The  diagonal  elements  of  the  covariance  matrix  obtained  in  the  previous  section  contain  the 
variance  of  the  residual  stress  at  each  depth.  Approximate  95%  lower  and  upper  bound  of  residual 
stress  at  each  depth  can  be  obtained  from  Equation  38 

Upper  Bound  —  RSj  +  1.96(7/ 

Lower  Bound  —  RSj  —  1.96(7,  (38) 

where  RSj  is  the  residual  stress  obtained  from  the  deterministic  analysis  at  a  depth  i  and  o,  is  the 
standard  deviation  at  a  depth  i.  The  index  i  ranges  from  0  (surface)  to  a  depth  of  interest.  The 
number  1.96  represents  the  95%  bounds  for  normal  distribution.  The  upper  and  lower  confidence 
band  for  the  entire  depth  is  obtained  by  interpolating  the  upper  and  lower  confidence  bounds 
respectively.  The  confidence  band  implies  that  the  true  residual  stress  field  will  be  between  the 
confidence  band  95%  of  the  times. 

Validation  of  Confidence  Band 

A  Monte  Carlo  analysis  is  performed  to  demonstrate  the  validation  of  the  confidence  band 
developed.  Figure  64  shows  the  flow  chart  of  the  validation  procedure.  Random  samples  of  error 
are  generated  from  a  normal  distribution  with  a  zero  mean  and  a  variance  of  mean  squared  error 
obtained  from  the  non-linear  regression  analysis.  The  deterministic  fit  is  added  to  the  random 
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Figure  65:  JC  Model  Fit  to  the  Experimental  Data 

error.  New  estimates  of  model  constants  are  obtained  from  non-linear  regression  analysis.  The 
developed  methodology  is  used  to  obtain  the  confidence  band  for  the  entire  depth  of  interest  for 
the  new  simulated  experimental  data.  The  deterministic  residual  stress  is  verified  at  each  depth 
and  if  the  stress  lies  in  between  the  confidence  bounds,  a  counter  is  prepared.  This  procedure  is 
repeated  several  times  and  the  counter  is  updated  each  time  the  deterministic  residual  stress  is  in 
between  the  confidence  bounds.  The  counter  at  each  depth  can  be  used  to  approximate  the  true 
confidence  level.  A  plot  of  obtained  confidence  level  vs.  depth  is  prepared  to  compare  with  the 
target  level. 

6.4  Regression  (Input)  Uncertainty  Quantification  Results 

Experimental  data  [75]  of  material  behavior  for  Ti-6A1-4V  at  four  different  strain  rates  are  shown 
in  Figure  65.  This  figure  shows  the  plot  of  flow  stress  vs.  plastic  strain.  The  experiments  have 
been  performed  using  the  Hopkinson  bar  test,  details  can  be  found  in  the  cited  reference.  The  JC 
model  is  used  to  curve  fit  the  available  experimental  data  and  will  later  be  used  to  predict  the 
material  behavior  at  higher  strain  rates.  The  JC  model  details  are  provided  in  Chapter  3  and  is 
shown  in  Equation  39  for  continuity.  Figure  65  shows  that  the  JC  model  was  able  to  match  the 
experimental  data  at  three  strain  rates  of  20.v  1 ,  1.5s-1,  and  0.045“ 1  relatively  well  compared  to 
strain  rate  of  1.5s-1.  Model  constants  A,  B ,  n,  and  C  are  estimated  using  non-linear  regression 
analysis.  Co  is  assumed  to  be  equal  to  Is-1. 

e  =  {A+Ben]  1  +C/n—  (39) 

L  £o_ 

Figure  66  shows  another  representation  of  the  fit  to  the  experimental  data.  Comparison  between 
actual  values  vs.  predicted  values  are  plotted  at  different  actual  values.  A  plot  of  straight  line  with 
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Figure  66:  JC  Model  Fit  vs  Actual  Experimental  Data 

45°  angle  represents  a  good  fit  with  the  experimental  data  at  all  predicted  values.  Figure  66  shows 
a  straight  line  with  better  fit  to  the  experimental  data  at  lower  predicted  values  compared  to 
higher.  Two  assumptions  of  random  errors:  (i)  constant  variance  and  (ii)  zero  mean  can  be 
checked  at  the  same  time  using  the  residual  plot.  A  residual  plot  is  a  plot  of  residuals  against  the 
fitted  values.  If  both  the  assumptions  are  satisfied  the  scatter  is  expected  to  vary  randomly  around 
zero  for  all  predicted  values.  Figure  67  shows  the  residual  plot  for  the  JC  model  from  the 
non-linear  regression  analysis.  The  figure  shows  that  both  the  assumptions  were  satisfied  for  all 
predicted  values. 

The  model  constant  estimates  obtained  are  shown  in  Table  24. 


Table  24:  Material  Model  Constant  Estimates 


Model  Constant  Estimates 

Value 

A 

1170.07  MPa 

B 

837.34  MPa 

n 

0.3409 

C 

0.017 

An  estimate  of  the  covariance  matrix  for  the  model  constant  estimates  is  obtained  from  the 
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Figure  67:  JC  Model  Fit  vs  Residuals 


non-linear  regression  analysis  using  Equation  20  and  is  shown  in  Equation  40: 

_  A  A  /V  _ 

A  B  n  C 

A  5.72  x  102  1.51  xlO3  1.18  -8.82  x  10~4 

B  7.10  xlO3  4.08  -4.83  x  10"3  (40) 

n  2.75  x  10"3  -3.58  x  10”6 

_  C  Sym  4.29  x  10^7 

The  variance  of  the  model  constant  estimates  are  along  the  diagonals.  The  off-diagonal  terms 
represent  the  covariance  between  model  constant  estimates.  A  correlation  coefficient  matrix  is 
obtained  to  provide  a  better  interpretation  of  the  covariance  matrix  because  the  model  constants 
are  on  different  scales.  Equation  41  shows  the  formula  to  obtain  the  correlation  coefficient 
between  two  random  variables,  X  and  Y : 


corr  [X  .Y] 


cov  [X.  Y] 


(41) 


where  cov  [X.Y]  is  the  covariance  between  X  and  Y  and  Var  [•]  represents  variance  of  a  variable. 
Equation  42  shows  the  calculated  correlation  coefficient  matrix: 

A  A  A  - 

A  B  n  C 

A  1  0.75  0.94  -0.06 

[corr]  -  B  1  0.92  -0.09  (42) 

n  1  -0.1 

_  C  Sym  1 
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Figure  68:  Bootstrap  Histogram  of  A 


Equation  42  indicates  that  the  estimate  of  model  constant  C  is  negatively  correlated  with  other 
estimates  and  also  these  correlations  are  low.  Model  constant  estimate  h  is  highly  correlated  with 

/v  /V 

estimates  A  and  B. 

Fixed-x  bootstrap  sampling  technique  is  performed  to  evaluate  the  multivariate  normality 
assumption  of  the  model  constant  estimates.  Figures  68,  69,  70,  and  71  show  the  bootstrap 
histograms  of  the  estimates  A,  B,  C,  and  h.  5000  bootstrap  estimates  are  obtained.  It  can  be 
seen  from  these  figures  that  the  model  constant  estimates  n  and  C  are  approximately  normal  while 

/•V  /V 

the  estimates  A  and  B  are  skewed.  The  bootstrap  method  provides  a  qualitative  validation  on  the 
multi-variate  normality  assumption  that  the  model  constant  estimates  closely  resemble  the 
assumed  distribution. 

6.5  Demonstration  Examples  and  LP  Application 

The  input  uncertainty  quantified  in  the  previous  section  is  propagated  through  closed-form 
demonstration  examples  and  the  FP  application  to  obtain  the  confidence  band.  Two  demonstration 
examples  are  investigated.  For  these  demonstration  examples,  an  explicit  relationship  between  the 
residual  stress  and  the  model  constant  estimate  is  assumed  such  that  the  residual  stress  field 
obtained  is  similar  to  the  residual  stress  field  generated  by  the  FP  application.  Monte  Carlo 
analysis  is  performed  for  the  demonstration  examples  to  validate  the  confidence  band. 

6.5.1  Demonstration  Example  1 

The  relationship  between  the  residual  stress  field  and  the  model  constant  estimates  is  shown  in 
Equation  43. 
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Figure  69:  Bootstrap  Histogram  of  B 
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Figure  70:  Bootstrap  Histogram  of  h 
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Figure  71:  Bootstrap  Histogram  of  C 


RS  =  — 1000  +  1000 log  (A  +  20 Bd  +  nd2  +  Cd 3)  (43) 

where  RS  is  the  residual  stress  value,  d  is  the  depth  in  mm,  and  A,  B,n,  and  C  are  the  model 
constants.  Equation  43  indicates  a  logarithmic  relationship  between  the  residual  stress  field  at  any 
depth  d  and  the  model  constants.  The  estimates  obtained  in  Table  24  are  used  to  obtain  the 
deterministic  residual  stress  field.  A  depth  of  2  mm  is  considered  and  80  equal  divisions  are 
assumed  representing  the  discretization  (FEA  mesh)  of  the  depth.  Figure  72  shows  the 
deterministic  residual  stress  field.  As  we  can  see  from  the  figure,  the  residual  stress  field  is 
compressive  at  the  surface  and  progressively  leads  to  tensile  at  sub  surface  which  is  similar  to  the 
laser  peened  residual  stress  field.  The  closed  form  equation  is  a  representative  of  the  residual 
stress  field  obtained  from  the  LP  simulation. 

A  central  difference  scheme  is  used  to  obtain  the  numerical  gradient  of  the  residual  stress  field 
with  respect  to  each  model  constant  estimate.  A  comparison  is  made  with  the  analytical  gradients 
because  the  explicit  relationship  between  the  residual  stress  field  and  the  model  constants  is 
known.  Figures  73,  74,  75,  and  76  show  the  comparison  between  analytical  and  numerical 
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gradient  of  residual  stress  field  for  model  constant  estimates  A,  B,  h,  and  C  respectively. 
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Deterministic  Residual  Stress  Field 


Figure  72:  Deterministic  Residual  Stress  Field  for  Demonstration  Problem  1 


Figure  73:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  A  for 
Demonstration  Problem  1 
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Figure  74:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  B  for 
Demonstration  Problem  1 


Figure  75:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  h  for 
Demonstration  Problem  1 
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Figure  76:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  C  for 
Demonstration  Problem  1 

A  good  agreement  between  the  numerical  and  analytical  gradient  is  evident  from  the  figures. 
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Increasing  trend  can  be  seen  for  the  model  constant  estimates  B,  C,  and  h,  while  decreasing  trend 
is  evident  for  A.  The  covariance  matrix  of  the  residual  stress  is  obtained  from  the  procedure 
described  in  Section  .  Figure  77  shows  the  95%  confidence  band  on  the  residual  stress  field  along 
with  the  deterministic  residual  stress  field. 
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Figure  77:  95%  Confidence  Band  on  the  Residual  Stress  Field  for  Demonstration  Problem  1 

The  confidence  band  indicates  that  the  true  residual  stress  field  is  in  between  the  band  95%  of  the 
times.  A  plot  of  standard  deviation  of  the  residual  stress  along  the  depth  is  shown  in  Figure  78. 


Figure  78:  Standard  Deviation  of  Residual  Stress  Field  for  Demonstration  Problem  1 
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A  lower  value  of  standard  deviation  can  be  seen  at  the  surface  and  progressively  increases  along 
the  depth.  The  smallest  variability  of  2  MPa  is  observed  at  the  surface  while  a  maximum 
variability  of  10  MPa  can  be  seen  at  a  depth  of  2  mm. 

A  MCS  on  Equation  43  is  performed  to  validate  the  developed  residual  stress  band.  Figure  79 
shows  the  comparison  between  the  target  and  achieved  confidence  level  for  the  million  Monte 
Carlo  samples  along  the  entire  depth. 


Figure  79:  Confidence  Level  Comparison  to  Target  Confidence  Level  for  Demonstration 
Problem  1 

A  lowest  confidence  level  of  «  92%  is  achieved  at  the  sub  surface  and  a  highest  confidence  of 
~  95%  at  near  surface.  Figure  79  indicates  that  the  developed  confidence  bands  are  able  to 
capture  the  target  confidence  levels  with  reasonable  accuracy. 

6.5.2  Demonstration  Example  2 

A  response  surface  fit  is  chosen  as  the  second  demonstration  example.  The  response  surface  is 
shown  in  Equation  44 


RS=  a  i  +  cc2Ada 3  +  aABdas  +  a6n  dai  +  cc^Cd**9  (44) 

where  i  =  1, 2, . .  .9  need  to  be  determined.  The  second  demonstration  example  represents  a 
polynomial  relationship  between  the  residual  stress  field  and  the  model  constant  estimates  as 
opposed  to  logarithmic  relationship  in  the  demonstration  example  1.  Also,  the  response  surface  is 
a  representative  of  the  residual  stress  field  induced  by  the  FP  process.  21  latin  hypercube  samples 
are  generated  using  JUMP  software  and  the  design  matrix  is  shown  in  Table  25.  Upper  and  lower 
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bounds  for  each  model  constant  estimate  is  provided.  The  design  matrix  is  generated  such  that  the 
samples  are  generated  from  the  entire  design  space.  FEA  is  performed  at  all  the  21  design  points 
and  the  residual  stress  values  are  obtained.  The  first  20  samples  are  used  to  obtain  the  coefficients 
and  the  21st  design  point  is  used  to  compare  the  response  surface  fit  and  the  FEA  result. 


Table  25:  Latin  Hypercube  Samples  Design  Matrix 


# 

A 

B 

n 

C 

d 

1 

1168.3 

862.59 

0.3881 

0.0174 

0.3 

2 

1180.2 

921.58 

0.3723 

0.0178 

1.4 

3 

1189.8 

803.6 

0.3776 

0.0168 

0.6 

4 

1173.1 

837.31 

0.3933 

0.0173 

1.7 

5 

1185.0 

904.73 

0.2989 

0.0177 

0.8 

6 

1177.8 

786.74 

0.3199 

0.0178 

1.1 

7 

1151.5 

761.46 

0.3304 

0.0172 

1.3 

8 

1175.4 

795.17 

0.2884 

0.0171 

1.6 

9 

1187.4 

913.15 

0.3566 

0.0169 

1.5 

10 

1161.1 

769.89 

0.3618 

0.0169 

0.2 

11 

1163.4 

845.73 

0.3356 

0.0167 

2.0 

12 

1170.7 

753.03 

0.3828 

0.0177 

0.5 

13 

1156.3 

871.02 

0.2936 

0.0174 

1.9 

14 

1192.2 

879.44 

0.3409 

0.0171 

0.0 

15 

1165.9 

887.87 

0.3146 

0.017 

0.7 

16 

1153.9 

812.02 

0.3461 

0.0179 

1.8 

17 

1158.7 

820.45 

0.3251 

0.0176 

0.1 

18 

1194.6 

854.16 

0.3041 

0.0168 

1.0 

19 

1182.6 

778.31 

0.3094 

0.0172 

0.4 

20 

1149.1 

896.3 

0.3513 

0.0175 

1.2 

21 

1197.0 

828.88 

0.3671 

0.0175 

0.9 

Figure  80  shows  the  comparison  between  the  FEA  result  and  response  surface  at  the  2 1 st  design 
point  up  to  a  depth  of  2  mm. 
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Figure  80:  Response  Surface  Fit  Comparison  with  FEA 


It  is  evident  from  Figure  80  that  the  response  surface  is  in  good  agreement  with  the  FEA  result 
along  the  entire  depth  except  at  the  surface.  The  response  surface  over  predicts  at  the  surface,  but 
for  the  purpose  of  validation  of  the  framework  the  above  obtained  response  surface  is  sufficient. 
The  coefficients  obtained  from  the  response  surface  fit  are  shown  in  Table  26. 


Table  26:  Response  Surface  Coefficients 


Parameter 

Value 

a\ 

-54.70 

a2 

-0.02 

«3 

3.45 

«4 

0.2304 

«5 

0.6984 

a6 

-345.64 

a2 

-0.5499 

«8 

155.3663 

CLg 

5.6861 

A  depth  of  2  mm  is  considered  and  80  equal  divisions  are  assumed  representing  the  discretization 
(FEA  mesh)  of  the  depth.  Figure  81  shows  the  deterministic  residual  stress  field. 
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Deterministic  Residual  Stress  Field 


Figure  81:  Deterministic  Residual  Stress  Field  for  Demonstration  Example  2 


As  we  can  see  from  the  figure,  the  same  trend  observed  in  demonstration  example  1  is  also  seen 
here.  The  residual  stress  field  consists  of  compressive  stresses  at  the  surface  and  progressively 
leads  to  tensile  stresses  at  subsurface  which  is  similar  to  the  LP  residual  stress  field.  A  higher 
compressive  residual  stress  of  1000  MPa  is  obtained  at  the  surface  compared  to  the  previous 
example  which  had  a  compressive  residual  stress  of  250  MPa  at  the  surface.  A  central  difference 
scheme  is  used  again  to  obtain  the  numerical  gradients  of  the  residual  stress  field  with  respect  to 
each  model  constant  estimate.  A  comparison  is  made  with  the  analytical  gradients  because  the 
explicit  relationship  between  the  residual  stress  field  and  the  model  constants  is  known.  Figures 
82,  83,  84,  and  85  show  the  comparison  between  analytical  and  numerical  gradient  of  residual 
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stress  field  for  model  constant  estimates  A,  B ,  n,  and  C  respectively. 
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Figure  82:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  A  for 
Demonstration  Problem  2 


Figure  83:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  B  for 
Demonstration  Problem  2 
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Figure  84:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  n  for 
Demonstration  Problem  2 
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Figure  85:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  C  for 
Demonstration  Problem  2 

A  good  agreement  between  the  numerical  and  analytical  gradient  is  evident  from  these  figures. 
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Increasing  trend  can  be  seen  for  the  model  constant  estimates  B,  C,  and  h,  while  decreasing  trend 
is  evident  for  A.  The  covariance  matrix  of  the  residual  stress  is  obtained  from  the  sensitivity 
matrix  of  the  residual  stress  field  with  respect  to  each  model  constant  estimate.  Figure  86  shows 
the  95%  confidence  band  on  the  residual  stress  field  along  with  the  deterministic  residual  stress 
field. 
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Figure  86:  95%  Confidence  Band  on  the  Residual  Stress  Field  for  Demonstration  Problem  2 


The  confidence  band  indicates  that  the  true  residual  stress  field  is  in  between  the  band  95%  of  the 
times.  A  plot  of  the  variability  of  the  residual  stress  field  along  the  depth  can  be  see  in  Figure  87 
that  show  standard  deviation  against  depth. 


Figure  87:  Standard  Deviation  Variation  along  the  Depth  for  Demonstration  Problem  2 
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A  higher  variability  of  residual  stress  can  be  seen  at  the  surface  and  progressively  decreases  along 
the  depth.  A  maximum  variability  of  ~  160MPa  is  obtained  at  the  surface  while  least  variability 
of  sM5MPa  can  be  seen  at  a  depth  of  ~  1  mm. 

A  Monte  Carlo  analysis  is  performed  to  validate  the  developed  residual  stress  band.  Figure  88 
shows  the  comparison  between  the  target  and  achieved  confidence  level  for  the  million  Monte 
Carlo  samples  along  the  entire  depth.  It  can  be  seen  from  Figure  88  that  the  achieved  confidence 
level  is  lower  than  the  target  confidence  level  at  the  surface  regions  while  a  higher  confidence 
level  is  achieved  sub  surface. 


Figure  88:  Confidence  Level  Comparison  to  Target  Confidence  Level  for  Demonstration 
Problem  2 

A  lowest  confidence  level  of  «  93%  is  achieved  at  a  depth  of  ~  0.6 mm  and  a  highest  confidence 
of  ~  96.5%  at  a  depth  of  ~  1  Amin.  The  obtained  confidence  level  differs  from  the  target  value  of 
95%  could  be  due  to  Taylor  series  approximation  and/or  lack  of  normality  in  regression  estimates. 

6.5.3  LP  Application 

For  the  LP  application,  no  explicit  relationship  between  residual  stress  and  material  model 
constants  is  available.  FEA  is  used  to  simulate  the  residual  stress  fields  generated  by  the  LP 
process.  The  estimates  obtained  in  Table  24  are  used  to  obtain  the  deterministic  residual  stress 
field.  Figure  89  shows  the  deterministic  residual  stress  field  with  a  peak  pressure  of  5.5  GPa. 
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Figure  89:  Deterministic  Residual  Stress  Field  for  the  LP  Process  Simulation 


A  comparison  is  also  shown  with  the  experimental  data  [38].  It  is  evident  from  the  Figure  89  that 
compressive  residual  stresses  are  generated  at  the  surface  and  compensated  by  tensile  stresses  in 
the  sub-surface.  Figures  90,  91,  92,  and  93  show  the  gradients  of  residual  stress  field  with  respect 
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to  model  constant  estimates  A,  B ,  h.  and  C  and  respectively. 
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Figure  90:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  A  for 
the  LP  Simulation 


Figure  91:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  B  for 
the  LP  Simulation 
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Figure  92:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  h  for 
the  LP  Simulation 
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9000 


Figure  93:  Gradient  of  Residual  Stress  Field  with  respect  to  Model  Constant  Estimate  C  for 
the  LP  Simulation 

Central  difference  scheme  is  used  to  obtain  the  gradients.  Assuming  the  true  gradients  are  smooth 
function,  the  numerical  gradients  generated  by  central  difference  scheme  are  smoothed.  It  can  be 
seen  from  these  figures  that  the  gradients  are  on  different  scales.  Decreasing  trend  can  be  seen  for 

yv  /v  /\ 

the  model  constant  estimates  A,  B,  and  C,  while  an  increasing  trend  is  evident  for  estimate  h.  95% 
confidence  band  for  the  residual  stress  field  at  each  depth  is  shown  in  Figure  94. 
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Figure  94:  95%  Confidence  Band  on  the  Residual  Stress  Field  for  the  LP  Simulation 

A  higher  variance  of  residual  stress  is  obtained  at  the  surface  and  decreases  along  the  depth.  A 
plot  of  the  variability  along  the  depth  can  be  seen  in  Figure  95  that  shows  the  standard  deviation 
against  depth. 


Figure  95:  Standard  Deviation  of  Residual  Stress  Field  for  the  LP  Simulation 
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A  maximum  deviation  of  7MPa  is  obtained  at  the  surface.  This  indicates  that  the  JC  model  is 
robust  in  prediction  of  residual  stress  field  for  the  LP  process.  The  variation  of  residual  stress 
depends  on  the  material  model,  experimental  data,  and  relationship  between  residual  stress  field 
and  model  constants. 

6.6  Section  Summary 

In  this  chapter,  a  framework  is  developed  to  quantify  the  uncertainty  in  residual  stress  field 
induced  by  the  LP  simulation.  Input  uncertainty  is  quantified  from  a  non-linear  regression 
analysis  that  is  used  to  fit  the  constitutive  material  model  to  the  experimental  data  of  stress-strain 
at  different  strain  rates.  A  bootstrapping  technique  is  used  to  evaluate  the  multivariate  normality 
assumption  of  the  material  model  constant  estimates.  Taylor  series  expansion  in  combination  with 
sensitivity  analysis  is  used  to  propagate  the  input  uncertainty.  Two  demonstration  examples  are 
shown  to  validate  the  methodology  by  comparing  the  obtained  confidence  level  with  target 
confidence  level  using  Monte  Carlo  analysis.  The  JC  material  model  is  shown  to  be  a  robust 
model  to  simulate  the  LP  process  with  a  maximum  residual  stress  standard  deviation  of  7  MPa  is 
obtained  at  the  surface. 
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7.0  SUMMARY  AND  FUTURE  DIRECTIONS 

This  chapter  provides  a  brief  summary  of  the  research  work  performed  and  possible  future 
directions  for  the  research  that  can  be  extended.  In  this  work,  the  simulation  of  the  LP  process  is 
performed  and  the  residual  stress  field  generated  by  the  LP  simulation  is  quantified  by  the 
propagation  of  regression  uncertainty. 

Material  model  validation  for  the  simulation  of  the  LP  process  is  presented.  FEA  is  shown  to  be  a 
useful  method  to  simulate  the  LP  process.  Three  different  material  models,  the  EPP,  the  JC 
model,  and  the  ZA  model,  are  used  for  the  validation  purposes.  Efficient  modeling  of  the  LP 
process  is  performed  by  incorporating  infinite  elements  to  represent  the  non-reflecting  boundary 
conditions.  An  integration  framework  is  implemented  for  efficient  data  management  between 
ABAQUS/Explicit  and  ABAQUS/Standard.  The  EPP  model,  which  is  most  often  used  in  the 
literature,  is  shown  to  produce  inconsistent  results.  The  simulation  results  are  compared  with  the 
experimental  data  for  three  different  peak  pressures  including  multiple  shots.  The  ZA  model, 
which  is  based  on  dislocation  mechanics,  produces  consistent  trends  but  overestimates  the  results 
compared  to  experimental  data.  The  JC  model  is  shown  to  produce  consistent  results  matching 
the  trends  and  better  agreement  with  experimental  results. 

Inverse  optimization-based  approach  is  designed  to  obtain  the  material  behavior  when  very  little 
or  no  experimental  data  of  stress-strain  curves  is  available.  The  optimization-based  approach  is 
shown  to  predict  residual  stresses  that  are  consistent  with  experimental  results.  The  consistency 
of  the  approach  is  shown  by  validating  for  two  materials  including  Inconel®718  and  Ti-6A1-4V. 
LP  experiments  were  performed  in  collaboration  with  LSP  Technologies  Inc,  Dublin,  OH  with  a 
Nd-glass  laser  for  Inconel®718  at  four  different  energy  densities  and  the  residual  stress 
measurements  were  made  using  an  x-ray  diffraction  method.  For  the  Inconel®718,  the  JC  and  the 
KHL  models  predicted  the  trends  and  the  simulation  results  are  in  agreement  with  the 
experimental  results  for  the  lower  two  peak  pressures.  The  JC  and  the  KHL  models  are  shown  to 
perform  better  than  the  ZA  model  in  prediction  of  residual  stresses  compared  to  experimental  data 
for  Ti-6A1-4V. 

A  framework  is  developed  to  quantify  the  uncertainty  in  residual  stress  field  induced  by  the  LP 
simulation.  The  input  uncertainty  is  quantified  from  a  non-linear  regression  analysis  that  is  used 
to  fit  the  constitutive  material  model  to  the  experimental  data  of  stress-strain  at  different  strain 
rates.  A  technique  from  statistics  known  as  ’bootstrap  for  regression’  is  used  to  evaluate  the 
multivariate  normality  assumption  of  the  material  model  constant  estimates.  Taylor  series 
expansion  in  combination  with  sensitivity  analysis  is  used  to  propagate  the  input  uncertainty.  Two 
demonstration  examples  are  shown  to  validate  the  methodology  by  comparing  the  obtained 
confidence  level  with  target  confidence  level  using  Monte  Carlo  analysis.  The  first  example  is  a 
logarithmic  relationship  between  the  residual  stress  field  and  the  model  constant  estimates  while 
the  second  example  is  a  response  surface  fit  of  the  FEA.  The  LP  application  results  show  that  the 
JC  material  model  is  a  robust  model  to  simulate  the  residual  stresses  induced  by  the  LP  process. 
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7.1  Future  Directions 


The  presented  research  work  of  finite  element  simulations  of  the  residual  stresses  induced  by  the 
LP  process  can  be  extended  to  modeling  of  thin  specimens,  two  side  LP  simulation,  geometric 
effects  of  LP,  and  integrated  framework  for  LP  design. 

7.1.1  Modeling  of  Thin  Specimens 

In  the  current  work,  LP  simulations  are  performed  on  thick  specimens.  Future  work  can  be 
extending  the  research  work  to  thin  specimens.  In  this  aspect,  accurate  tracking  of  shock  waves 
reflecting  from  the  surface  and  appropriate  material  behavior  changes  can  be  investigated. 

7.1.2  Two  Sided  LP  Simulation 

Very  little  work  has  been  performed  in  the  literature  on  two  sided  LP  simulation.  Simulating  the 
two  sided  LP  process  can  be  a  challenge  not  only  in  understanding  the  physics  of  the  material  but 
also  in  computational  cost.  Efficient  techniques  need  to  be  developed  to  address  the  computing 
cost. 

7.1.3  Geometric  Effects  of  the  LP  Process 

In  the  current  work  and  most  of  the  literature  focuses  the  work  on  simulating  flat  surfaces. 
Geometric  effects  of  component  that  need  to  be  laser  peened  can  be  investigated.  Efforts  must  be 
made  to  model  the  real  scale  component  instead  of  coupon  level  scale,  that  will  allow  curvature 
effects  to  be  investigated. 

7.1.4  Bootstrapping  Method  for  Parametric  Uncertainty 

In  this  work,  the  bootstrapping  method  was  used  to  validate  normality  assumption  of  the  model 
constant  estimates  for  the  non-linear  regression  analysis.  For  the  future  work,  the  bootstrapping 
method  can  be  adapted  to  define  the  distribution  of  a  parametric  uncertainty  including  peak 
pressure,  spot  radius,  material  properties,  and  geometric  dimensions  .  Figure  96  shows  the 
schematic  representation  of  application  of  bootstrapping  technique  for  regression  uncertainty. 
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Figure  96:  Bootstrapping  for  Parametric  Uncertainty 

The  peak  pressure  of  the  shock  wave  can  be  considered  as  an  example  for  demonstration.  The 
bootstrapping  method  can  be  used  to  obtain  a  histogram  of  the  peak  pressure  with  limited 
experimental  data.  The  mean  of  the  experimental  data  is  considered  as  the  sample  statistic  to 
obtain  the  histogram.  The  obtained  histogram  is  matched  with  the  known  distributions  to  obtain 
not  only  the  nature  of  the  distribution  but  also  the  parameters  of  the  distribution.  Other  sample 
statistics  such  as  variance  could  also  be  used  to  obtain  further  analysis  of  the  parameters  of  the 
distribution. 

Other  types  of  uncertainties  including  model  uncertainty  and  shape  uncertainty  can  be  integrated 
with  the  bootstrapping  method  based  on  the  availability  of  the  experimental  data.  For  the  model 
uncertainty,  bootstrapping  method  can  be  used  to  obtain  the  confidence  interval  of  the  LP  residual 
stresses  and  then  integrated  with  the  model  uncertainty  quantification  method.  Bootstrapping 
method  can  be  used  with  multiple  statistics  including  mean  and  standard  deviation  to  define  the 
shape  uncertainty. 

7.1.5  Integration  Framework  for  LP  Design 

An  integrated  framework  can  be  investigated  for  the  implementation  of  LP  design  as  shown  in  the 
Figure  97.  Different  tools  can  be  integrated  to  obtain  maximum  benefit  from  the  LP  process. 
Given  a  structural  component,  this  framework  would  consider  both  its  geometric  configuration 
and  constituent  material(s)  to  determine  an  optimal  LP  configuration.  The  existing  database  of 
experimental  work  based  on  geometry,  LP  parameters,  and  materials  would  be  used  to  obtain  LP 
process  parameters.  High  fidelity  simulations  combined  with  semi-empirical  relations  have  to  be 
developed  and  validated  for  new  structures  to  update  this  database.  The  fatigue  life  estimation 
methodology  including  stress  relaxation  mechanisms  have  to  be  implemented  simultaneously  to 
obtain  a  realistic  life  estimate  of  the  structural  component.  Vibration  characteristics  can  also  be 
investigated  to  avoid  resonant  frequencies.  Risk-based  analysis  and  optimization  of  the  LP 
process  would  be  performed  based  on  experimental  and  simulated  data  to  obtain  risk  quantified 
optimum  design.  This  procedure  would  then  be  iterated  for  a  specific  design  until  a  user  input 
criterion  was  met.  Examples  include  a  percentage  increase  in  fatigue  life  or  failure  at  a  given 
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location.  Eight  disciplines  are  represented  in  the  figure,  while  other  disciplines  can  be  added  to 
the  framework. 
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Figure  97:  Integration  Framework  for  LP  Design 
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